From 0b892dbe3155711c4560e61832f10f4a0fcf768a Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 25 Oct 2023 21:38:20 +0200 Subject: [PATCH] restructuration QuAcK for restricted and unrestricted branches --- input/methods | 6 +- input/options | 24 +- src/CC/{CC.f90 => RCC.f90} | 4 +- src/CI/CID.f90 | 44 +-- src/CI/CISD.f90 | 56 ++-- src/CI/RCI.f90 | 101 ++++++ src/CI/{CI.f90 => UCI.f90} | 21 +- src/GF/RGF.f90 | 148 +++++++++ src/GF/{GF.f90 => UGF.f90} | 59 +--- src/GT/RGT.f90 | 174 +++++++++++ src/GT/{GT.f90 => UGT.f90} | 72 +---- src/GW/RGW.f90 | 172 +++++++++++ src/GW/{GW.f90 => UGW.f90} | 59 +--- src/HF/HF.f90 | 150 --------- src/HF/RHF.f90 | 8 +- src/HF/ROHF.f90 | 10 +- src/HF/UHF.f90 | 4 +- src/HF/mix_guess.f90 | 9 +- src/MP/RMP.f90 | 62 ++++ src/MP/{MP.f90 => UMP.f90} | 20 +- src/QuAcK/QuAcK.f90 | 412 +++++-------------------- src/QuAcK/RQuAcK.f90 | 322 +++++++++++++++++++ src/QuAcK/UQuAcK.f90 | 338 ++++++++++++++++++++ src/QuAcK/read_methods.f90 | 132 ++++---- src/QuAcK/read_options.f90 | 96 ++---- src/RPA/RRPA.f90 | 103 +++++++ src/RPA/{RPA.f90 => URPA.f90} | 39 +-- src/utils/orthogonalization_matrix.f90 | 5 +- src/utils/read_molecule.f90 | 28 +- 29 files changed, 1752 insertions(+), 926 deletions(-) rename src/CC/{CC.f90 => RCC.f90} (96%) create mode 100644 src/CI/RCI.f90 rename src/CI/{CI.f90 => UCI.f90} (75%) create mode 100644 src/GF/RGF.f90 rename src/GF/{GF.f90 => UGF.f90} (63%) create mode 100644 src/GT/RGT.f90 rename src/GT/{GT.f90 => UGT.f90} (56%) create mode 100644 src/GW/RGW.f90 rename src/GW/{GW.f90 => UGW.f90} (63%) delete mode 100644 src/HF/HF.f90 create mode 100644 src/MP/RMP.f90 rename src/MP/{MP.f90 => UMP.f90} (73%) create mode 100644 src/QuAcK/RQuAcK.f90 create mode 100644 src/QuAcK/UQuAcK.f90 create mode 100644 src/RPA/RRPA.f90 rename src/RPA/{RPA.f90 => URPA.f90} (62%) diff --git a/input/methods b/input/methods index 5f36a66..b1f599f 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ -# RHF UHF ROHF RMOM UMOM KS - F T F F F F +# RHF UHF GHF ROHF + T F F F # MP2* MP3 F F # CCD pCCD DCD CCSD CCSD(T) @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # 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 F F F F F F # * unrestricted version available diff --git a/input/options b/input/options index e73a4a6..1f2fdc6 100644 --- a/input/options +++ b/input/options @@ -1,17 +1,17 @@ -# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability - 5000 0.0000001 T 5 1 1 T 0.0 T +# HF: maxSCF thresh DIIS guess mix_guess level_shift stability + 5000 0.0000001 5 1 0.0 0.0 T # MP: reg F -# CC: maxSCF thresh DIIS n_diis - 64 0.0000001 T 5 -# spin: TDA singlet triplet spin_conserved spin_flip - F T T T T -# GF: maxSCF thresh DIIS n_diis lin eta renorm reg - 256 0.00001 T 5 T 0.0 0 F -# GW: maxSCF thresh DIIS n_diis lin eta TDA_W reg - 256 0.00001 T 5 F 0.0 F F -# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 256 0.00001 T 5 F 0.0 F F +# CC: maxSCF thresh DIIS + 64 0.0000001 5 +# spin: TDA spin_conserved spin_flip + F T T +# GF: maxSCF thresh DIIS lin eta renorm reg + 256 0.00001 5 F 0.0 0 F +# GW: maxSCF thresh DIIS lin eta TDA_W reg + 256 0.00001 5 F 0.0 F F +# GT: maxSCF thresh DIIS lin eta TDA_T reg + 256 0.00001 5 F 0.0 F F # ACFDT: AC Kx XBS F F T # BSE: phBSE phBSE2 ppBSE dBSE dTDA diff --git a/src/CC/CC.f90 b/src/CC/RCC.f90 similarity index 96% rename from src/CC/CC.f90 rename to src/CC/RCC.f90 index ddff7e5..2e5722d 100644 --- a/src/CC/CC.f90 +++ b/src/CC/RCC.f90 @@ -1,5 +1,5 @@ -subroutine CC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) +subroutine RCC(doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & + maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) ! Coupled-cluster module diff --git a/src/CI/CID.f90 b/src/CI/CID.f90 index a640cb8..df14e0b 100644 --- a/src/CI/CID.f90 +++ b/src/CI/CID.f90 @@ -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 @@ -14,7 +14,7 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi integer,intent(in) :: nOin integer,intent(in) :: nVin 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) :: E0 @@ -26,7 +26,7 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi integer :: nV integer :: nR - double precision,allocatable :: F(:,:) + double precision,allocatable :: eHF(:) double precision,allocatable :: sERI(:,:,:,:) double precision,allocatable :: ERI(:,:,:,:) @@ -61,9 +61,9 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi nV = 2*nVin nR = 2*nRin - allocate(F(nBas,nBas),sERI(nBas,nBas,nBas,nBas)) - - call spatial_to_spin_fock(nBasin,Fin,nBas,F) + allocate(eHF(nBas),sERI(nBas,nBas,nBas,nBas)) + + call spatial_to_spin_MO_energy(nBasin,eHFin,nBas,eHF) call spatial_to_spin_ERI(nBasin,ERIin,nBas,sERI) ! Antysymmetrize ERIs @@ -144,22 +144,22 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi kcld = kcld + 1 tmp = & 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) & - - F(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) & - + F(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) & - + F(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) & - - F(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) & - - F(nO+a,nO+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) & - + F(nO+a,nO+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) & - + F(nO+b,nO+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) & - - F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + - eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & + - eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & + + eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + - eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & + + eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + + eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + - eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & + + eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*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) & + - eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*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) & + - eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*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) & + + eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*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,c)*Kronecker_delta(b,d) & + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90 index 039d6b4..fdd928e 100644 --- a/src/CI/CISD.f90 +++ b/src/CI/CISD.f90 @@ -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 @@ -14,7 +14,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI integer,intent(in) :: nOin integer,intent(in) :: nVin 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) :: E0 @@ -26,7 +26,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI integer :: nV integer :: nR - double precision,allocatable :: F(:,:) + double precision,allocatable :: eHF(:) double precision,allocatable :: sERI(:,:,:,:) double precision,allocatable :: ERI(:,:,:,:) @@ -62,9 +62,9 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI nV = 2*nVin 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) ! Antysymmetrize ERIs @@ -111,7 +111,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI do a=1,nV-nR ia = ia + 1 - tmp = F(i,nO+a) + tmp = 0d0 H(ishift+1,jshift+ia) = 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 kc = kc + 1 - tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) & - - F(i,k)*Kronecker_delta(a,c) & - + F(nO+a,nO+c)*Kronecker_delta(i,k) & + tmp = E0*Kronecker_delta(i,k)*Kronecker_delta(a,c) & + - eHF(i)*Kronecker_Delta(i,k)*Kronecker_delta(a,c) & + + eHF(a)*Kronecker_delta(a,c)*Kronecker_delta(i,k) & - ERI(nO+a,k,nO+c,i) 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 kcld = kcld + 1 - tmp = - F(l,nO+d)*Kronecker_delta(a,c)*Kronecker_delta(i,k) & - + 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) & + tmp = - ERI(k,l,nO+d,i)*Kronecker_delta(a,c) & + 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,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 tmp = & 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) & - - F(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) & - + F(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) & - + F(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) & - - F(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) & - - F(nO+a,nO+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) & - + F(nO+a,nO+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) & - + F(nO+b,nO+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) & - - F(nO+b,nO+c)*Kronecker_delta(a,d)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & + - eHF(j)*Kronecker_delta(l,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,k) & + - eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(i,l) & + + eHF(j)*Kronecker_delta(k,j)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(i,l) & + - eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,k) & + + eHF(i)*Kronecker_delta(l,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,k) & + + eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*Kronecker_delta(j,l) & + - eHF(i)*Kronecker_delta(k,i)*Kronecker_delta(a,c)*Kronecker_delta(b,d)*Kronecker_delta(j,l) & + + eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*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) & + - eHF(a)*Kronecker_delta(a,d)*Kronecker_delta(b,c)*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) & + - eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*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) & + + eHF(b)*Kronecker_delta(b,d)*Kronecker_delta(a,c)*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,c)*Kronecker_delta(b,d) & + ERI(nO+a,l,nO+d,j)*Kronecker_delta(b,c)*Kronecker_delta(i,k) & diff --git a/src/CI/RCI.f90 b/src/CI/RCI.f90 new file mode 100644 index 0000000..e6e21f7 --- /dev/null +++ b/src/CI/RCI.f90 @@ -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 diff --git a/src/CI/CI.f90 b/src/CI/UCI.f90 similarity index 75% rename from src/CI/CI.f90 rename to src/CI/UCI.f90 index 0dfa91a..8ddae0e 100644 --- a/src/CI/CI.f90 +++ b/src/CI/UCI.f90 @@ -1,6 +1,5 @@ -subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin_conserved,spin_flip, & - nBas,nC,nO,nV,nR,nS,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int,dipole_int_aa,dipole_int_bb, & - epsHF,EHF,cHF,S,F) +subroutine 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,F) ! Configuration interaction module @@ -14,10 +13,7 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin logical :: doCID logical :: doCISD logical :: doFCI - logical :: unrestricted - logical,intent(in) :: singlet - logical,intent(in) :: triplet logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip 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) :: F(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_aabb(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_bb(nBas,nBas,ncart,nspin) @@ -50,12 +44,8 @@ subroutine CI(doCIS,doCIS_D,doCID,doCISD,doFCI,unrestricted,singlet,triplet,spin if(doCIS) then call wall_time(start_CI) - if(unrestricted) then - 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) - else - call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) - end if + 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) call wall_time(end_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 call wall_time(start_CI) - call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF) call wall_time(end_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 call wall_time(start_CI) - call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI,F,EHF) call wall_time(end_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) 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 diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 new file mode 100644 index 0000000..6ff348e --- /dev/null +++ b/src/GF/RGF.f90 @@ -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 diff --git a/src/GF/GF.f90 b/src/GF/UGF.f90 similarity index 63% rename from src/GF/GF.f90 rename to src/GF/UGF.f90 index 90174d2..1aa8357 100644 --- a/src/GF/GF.f90 +++ b/src/GF/UGF.f90 @@ -1,7 +1,7 @@ -subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,thresh,max_diis, & - dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,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, & - dipole_int_AO,dipole_int,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) +subroutine UGF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,renorm,maxSCF,thresh,max_diis, & + 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_aaaa,ERI_aabb,ERI_bbbb, & + dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) ! Green's function module @@ -15,7 +15,6 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t logical :: doqsGF2 logical :: doG0F3 logical :: doevGF3 - logical :: unrestricted integer :: renorm 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) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: singlet - logical,intent(in) :: triplet logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip 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) :: 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) :: ERI_aaaa(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) :: 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_bb(nBas,nBas,ncart) @@ -76,14 +71,9 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doG0F2) then call wall_time(start_GF) - if(unrestricted) then - 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, & - 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 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, & + dipole_int_aa,dipole_int_bb,epsHF) call wall_time(end_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 call wall_time(start_GF) - if(unrestricted) then - 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, & - 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 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, & + dipole_int_aa,dipole_int_bb,cHF,epsHF) call wall_time(end_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 call wall_time(start_GF) - if(unrestricted) then - 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, & - 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 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, & + 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 @@ -146,11 +125,7 @@ subroutine GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm,maxSCF,t if(doG0F3) then call wall_time(start_GF) - if(unrestricted) then - print*,'Unrestricted version of G0F3 not yet implemented! Sorry.' - else - call G0F3(renorm,nBas,nC,nO,nV,nR,ERI,epsHF) - end if + print*,'Unrestricted version of G0F3 not yet implemented! Sorry.' call wall_time(end_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 call wall_time(start_GF) - if(unrestricted) then - 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 + print*,'Unrestricted version of evGF3 not yet implemented! Sorry.' call wall_time(end_GF) t_GF = end_GF - start_GF diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 new file mode 100644 index 0000000..22eddf0 --- /dev/null +++ b/src/GT/RGT.f90 @@ -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 diff --git a/src/GT/GT.f90 b/src/GT/UGT.f90 similarity index 56% rename from src/GT/GT.f90 rename to src/GT/UGT.f90 index 74c765f..15c43ea 100644 --- a/src/GT/GT.f90 +++ b/src/GT/UGT.f90 @@ -1,8 +1,7 @@ -subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted,maxSCF,thresh,max_diis,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,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,dipole_int_AO,dipole_int,dipole_int_aa,dipole_int_bb, & - PHF,cHF,epsHF) +subroutine UGT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,maxSCF,thresh,max_diis,doACFDT, & + 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, & + ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) ! T-matrix module @@ -17,7 +16,6 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted logical :: doG0T0eh logical :: doevGTeh logical :: doqsGTeh - logical :: unrestricted integer,intent(in) :: maxSCF 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) :: dTDA logical,intent(in) :: doppBSE - logical,intent(in) :: singlet - logical,intent(in) :: triplet logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip 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) :: 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(inout):: ERI_aaaa(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(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_bb(nBas,nBas,ncart) @@ -83,14 +77,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(doG0T0pp) then call wall_time(start_GT) - if(unrestricted) then - 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, & - 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 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, & + dipole_int_aa,dipole_int_bb,cHF,epsHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -106,14 +95,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(doevGTpp) then 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, & - 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) - 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 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, & + cHF,epsHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -129,15 +113,9 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(doqsGTpp) then 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, & - 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) - 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 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, & + dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -153,14 +131,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(doG0T0eh) then call wall_time(start_GT) - if(unrestricted) then - 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 + print*,'Unrestricted version of G0T0eh not yet implemented! Sorry.' call wall_time(end_GT) t_GT = end_GT - start_GT @@ -176,12 +147,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(doevGTeh) then call wall_time(start_GT) - if(unrestricted) then - 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 + print*,'Unrestricted version of evGTeh not yet implemented! Sorry.' call wall_time(end_GT) t_GT = end_GT - start_GT @@ -197,13 +163,7 @@ subroutine GT(doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,unrestricted if(doqsGTeh) then call wall_time(start_GT) - if(unrestricted) then - 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 + print*,'Unrestricted version of qsGTeh not yet implemented! Sorry.' call wall_time(end_GT) t_GT = end_GT - start_GT diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 new file mode 100644 index 0000000..6a9ae9e --- /dev/null +++ b/src/GW/RGW.f90 @@ -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 diff --git a/src/GW/GW.f90 b/src/GW/UGW.f90 similarity index 63% rename from src/GW/GW.f90 rename to src/GW/UGW.f90 index 33696e2..70c810a 100644 --- a/src/GW/GW.f90 +++ b/src/GW/UGW.f90 @@ -1,8 +1,7 @@ -subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF,thresh,max_diis,doACFDT, & - exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet,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,dipole_int_AO,dipole_int,dipole_int_aa,dipole_int_bb, & - PHF,cHF,epsHF) +subroutine UGW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,maxSCF,thresh,max_diis,doACFDT, & + 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, & + ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) ! GW module @@ -17,7 +16,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF logical :: doufG0W0 logical :: doufGW logical :: doSRGqsGW - logical :: unrestricted integer,intent(in) :: maxSCF 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) :: dTDA logical,intent(in) :: doppBSE - logical,intent(in) :: singlet - logical,intent(in) :: triplet logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip 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) :: 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(inout):: ERI_aaaa(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(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_bb(nBas,nBas,ncart) @@ -82,14 +76,9 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF if(doG0W0) then call wall_time(start_GW) - if(unrestricted) then - 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, & - 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 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, & + dipole_int_aa,dipole_int_bb,cHF,epsHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -105,14 +94,9 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF if(doevGW) then call wall_time(start_GW) - if(unrestricted) then - 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, & - 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 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, & + ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,epsHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -128,16 +112,9 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF if(doqsGW) then 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, & - 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) - 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 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, & + dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -153,13 +130,7 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF if(doSRGqsGW) then call wall_time(start_GW) - if(unrestricted) then - 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 + print*,'Unrestricted version of SRG-qsGW not yet implemented! Sorry.' call wall_time(end_GW) t_GW = end_GW - start_GW @@ -175,7 +146,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF 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 @@ -191,7 +161,6 @@ subroutine GW(doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,unrestricted,maxSCF 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 diff --git a/src/HF/HF.f90 b/src/HF/HF.f90 deleted file mode 100644 index c26e0e1..0000000 --- a/src/HF/HF.f90 +++ /dev/null @@ -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 diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 034a383..5f61db6 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,5 +1,5 @@ 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 @@ -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 :: K(:,:) double precision,allocatable :: cp(:,:) + double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) ! 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) :: c(nBas,nBas) double precision,intent(out) :: P(nBas,nBas) - double precision,intent(out) :: F(nBas,nBas) ! Hello world @@ -73,8 +73,8 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc ! Memory allocation - allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),Fp(nBas,nBas), & - error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), & + Fp(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Guess coefficients and density matrix diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index 69a2701..ee24b32 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -1,5 +1,5 @@ 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 @@ -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) :: max_diis integer,intent(in) :: guess_type - logical,intent(in) :: mix + double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh 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 :: Fp(:,:) double precision,allocatable :: Ftot(:,:) - double precision,allocatable :: Ptot(:,:) + double precision,allocatable :: P(:,:,:) double precision,allocatable :: K(:,:,:) double precision,allocatable :: err(:,:) 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) :: e(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 @@ -79,7 +79,7 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc ! Memory allocation 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)) ! Guess coefficients and demsity matrices diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index a15dacb..87f37e2 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -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) :: max_diis integer,intent(in) :: guess_type - logical,intent(in) :: mix + double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh 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 - if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c) + if(nSCF == 1) call mix_guess(nBas,nO,mix,c) ! Compute density matrix diff --git a/src/HF/mix_guess.f90 b/src/HF/mix_guess.f90 index 9fc2ec0..5d037fe 100644 --- a/src/HF/mix_guess.f90 +++ b/src/HF/mix_guess.f90 @@ -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 @@ -9,15 +9,16 @@ subroutine mix_guess(nBas,nO,c) integer,intent(in) :: nBas integer,intent(in) :: nO(nspin) + double precision,intent(in) :: mix ! Local variables - double precision,parameter :: th = 0.1d0 - double precision,allocatable :: HOMO(:), LUMO(:) double precision,allocatable :: HOMOa(:),LUMOa(:) double precision,allocatable :: HOMOb(:),LUMOb(:) + double precision :: th + ! Output variables 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 + th = 2d0*pi*mix + HOMO(:) = c(:,nO(1),1) LUMO(:) = c(:,nO(1)+1,1) diff --git a/src/MP/RMP.f90 b/src/MP/RMP.f90 new file mode 100644 index 0000000..032ba9f --- /dev/null +++ b/src/MP/RMP.f90 @@ -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 diff --git a/src/MP/MP.f90 b/src/MP/UMP.f90 similarity index 73% rename from src/MP/MP.f90 rename to src/MP/UMP.f90 index c08a52b..14edae4 100644 --- a/src/MP/MP.f90 +++ b/src/MP/UMP.f90 @@ -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 @@ -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) :: doMP3 - logical,intent(in) :: unrestricted logical,intent(in) :: regularize 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) 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) + double precision,intent(in) :: epsHF(nBas,nspin) 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_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 call wall_time(start_MP) - if(unrestricted) then - 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 UMP2(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 @@ -58,13 +52,7 @@ subroutine MP(doMP2,doMP3,unrestricted,regularize,nBas,nC,nO,nV,nR,ERI,ERI_aaaa, if(doMP3) then call wall_time(start_MP) - - if(unrestricted) then - write(*,*) 'MP3 NYI for UHF reference' - stop - else - call MP3(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,epsHF) - end if + write(*,*) 'MP3 NYI for UHF reference' call wall_time(end_MP) t_MP = end_MP - start_MP diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 21ed7a1..e0801c1 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -3,35 +3,27 @@ program QuAcK implicit none include 'parameters.h' - logical :: unrestricted = .false. - logical :: doHF,doRHF,doUHF,doROHF,doRMOM,doUMOM + logical :: doRQuAcK,doUQuAcK,doGQuAcK + logical :: doRHF,doUHF,doGHF,doROHF logical :: dostab - logical :: doKS - logical :: doMP,doMP2,doMP3 - logical :: doCC,doCCD,dopCCD,doDCD,doCCSD,doCCSDT + logical :: doMP2,doMP3 + logical :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT logical :: dodrCCD,dorCCD,docrCCD,dolCCD - logical :: doCI,doCIS,doCIS_D,doCID,doCISD,doFCI - logical :: doRPA,dophRPA,dophRPAx,docrRPA,doppRPA - logical :: doGF,doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 - logical :: doGW,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW - logical :: doGT,doG0T0pp,doevGTpp,doqsGTpp - logical :: doG0T0eh,doevGTeh,doqsGTeh + logical :: doCIS,doCIS_D,doCID,doCISD,doFCI + logical :: dophRPA,dophRPAx,docrRPA,doppRPA + logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 + logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW + logical :: doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh integer :: nNuc,nBas - integer :: nEl(nspin) integer :: nC(nspin) integer :: nO(nspin) integer :: nV(nspin) integer :: nR(nspin) integer :: nS(nspin) - double precision :: ENuc,EHF + double precision :: ENuc 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 :: T(:,:) @@ -39,65 +31,41 @@ program QuAcK double precision,allocatable :: Hc(:,:) double precision,allocatable :: X(:,:) 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_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_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 + double precision :: start_QuAcK,end_QuAcK,t_QuAcK + double precision :: start_int ,end_int ,t_int integer :: maxSCF_HF,max_diis_HF - double precision :: thresh_HF,level_shift - logical :: DIIS_HF,guess_type,ortho_type,mix + double precision :: thresh_HF,level_shift,mix + logical :: guess_type logical :: reg_MP integer :: maxSCF_CC,max_diis_CC double precision :: thresh_CC - logical :: DIIS_CC - logical :: singlet - logical :: triplet logical :: spin_conserved logical :: spin_flip logical :: TDA - + integer :: maxSCF_GF,max_diis_GF,renorm_GF double precision :: thresh_GF - logical :: DIIS_GF,lin_GF,reg_GF + logical :: lin_GF,reg_GF double precision :: eta_GF integer :: maxSCF_GW,max_diis_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 integer :: maxSCF_GT,max_diis_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 logical :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA - + logical :: doACFDT,exchange_kernel,doXBS !-------------! ! Hello World ! @@ -123,31 +91,31 @@ program QuAcK ! Method selection ! !------------------! - call read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & - 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,doSRGqsGW, & - doufG0W0,doufGW, & - doG0T0pp,doevGTpp,doqsGTpp, & + call read_methods(doRHF,doUHF,doGHF,doROHF, & + 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,doSRGqsGW, & + doufG0W0,doufGW, & + doG0T0pp,doevGTpp,doqsGTpp, & doG0T0eh,doevGTeh,doqsGTeh) !--------------------------! ! Read options for methods ! !--------------------------! - call read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & - reg_MP, & - maxSCF_CC,thresh_CC,DIIS_CC,max_diis_CC, & - TDA,singlet,triplet,spin_conserved,spin_flip, & - maxSCF_GF,thresh_GF,DIIS_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_GT,thresh_GT,DIIS_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & - doACFDT,exchange_kernel,doXBS, & + call read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab, & + reg_MP, & + maxSCF_CC,thresh_CC,max_diis_CC, & + TDA,spin_conserved,spin_flip, & + maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, & + maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & + maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & + doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA) !------------------------------------------------! @@ -163,7 +131,7 @@ program QuAcK ! = nO*nV ! !------------------------------------------------! - call read_molecule(nNuc,nEl,nO,nC,nR) + call read_molecule(nNuc,nO,nC,nR) allocate(ZNuc(nNuc),rNuc(nNuc,ncart)) ! Read geometry @@ -183,9 +151,8 @@ program QuAcK ! 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), & - V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), & - dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas,nspin)) + allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), & + ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart)) ! Read integrals @@ -203,281 +170,56 @@ program QuAcK ! 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) - call HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,max_diis_HF, & - 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 + doGQuAcK = .false. + if(doGHF) doGQuAcK = .true. !-------------------------! -! 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) - call GF(doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,unrestricted,renorm_GF,maxSCF_GF,thresh_GF,max_diis_GF, & - dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,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_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_GF) + if(doUQuAcK) & + call 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) - t_GF = end_GF - start_GF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' - write(*,*) +!--------------------------! +! Generalized QuAcK branch ! +!--------------------------! - 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 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 +! if(doGQuAcK) call GQuAcK() !--------------! ! End of QuAcK ! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 new file mode 100644 index 0000000..796eb03 --- /dev/null +++ b/src/QuAcK/RQuAcK.f90 @@ -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 diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 new file mode 100644 index 0000000..e05de41 --- /dev/null +++ b/src/QuAcK/UQuAcK.f90 @@ -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 diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 0ccd26d..f20abcc 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,14 +1,14 @@ -subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & - doMP2,doMP3, & - doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & - do_drCCD,do_rCCD,do_crCCD,do_lCCD, & - doCIS,doCIS_D,doCID,doCISD,doFCI, & - dophRPA,dophRPAx,docrRPA,doppRPA, & - doG0F2,doevGF2,doqsGF2, & - doG0F3,doevGF3, & - doG0W0,doevGW,doqsGW,doSRGqsGW, & - doufG0W0,doufGW, & - doG0T0pp,doevGTpp,doqsGTpp, & +subroutine read_methods(doRHF,doUHF,doGHF,doROHF, & + doMP2,doMP3, & + doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & + do_drCCD,do_rCCD,do_crCCD,do_lCCD, & + doCIS,doCIS_D,doCID,doCISD,doFCI, & + dophRPA,dophRPAx,docrRPA,doppRPA, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & + doG0W0,doevGW,doqsGW,doSRGqsGW, & + doufG0W0,doufGW, & + doG0T0pp,doevGTpp,doqsGTpp, & doG0T0eh,doevGTeh,doqsGTeh) ! Read desired methods @@ -17,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & ! 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) :: doCCD,dopCCD,doDCD,doCCSD,doCCSDT 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 - character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6 + character(len=1) :: ans1,ans2,ans3,ans4,ans5,ans6 ! Open file with method specification @@ -40,19 +40,15 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doRHF = .false. doUHF = .false. + doGHF = .false. doROHF = .false. - doRMOM = .false. - doUMOM = .false. - doKS = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 - if(answer1 == 'T') doRHF = .true. - if(answer2 == 'T') doUHF = .true. - if(answer3 == 'T') doROHF = .true. - if(answer4 == 'T') doRMOM = .true. - if(answer5 == 'T') doUMOM = .true. - if(answer6 == 'T') doKS = .true. + read(1,*) ans1,ans2,ans3,ans4 + if(ans1 == 'T') doRHF = .true. + if(ans2 == 'T') doUHF = .true. + if(ans3 == 'T') doGHF = .true. + if(ans4 == 'T') doROHF = .true. ! Read MPn methods @@ -60,9 +56,9 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doMP3 = .false. read(1,*) - read(1,*) answer1,answer2 - if(answer1 == 'T') doMP2 = .true. - if(answer2 == 'T') doMP3 = .true. + read(1,*) ans1,ans2 + if(ans1 == 'T') doMP2 = .true. + if(ans2 == 'T') doMP3 = .true. ! Read CC methods @@ -73,12 +69,12 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doCCSDT = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 - if(answer1 == 'T') doCCD = .true. - if(answer2 == 'T') dopCCD = .true. - if(answer3 == 'T') doDCD = .true. - if(answer4 == 'T') doCCSD = .true. - if(answer5 == 'T') doCCSDT = .true. + read(1,*) ans1,ans2,ans3,ans4,ans5 + if(ans1 == 'T') doCCD = .true. + if(ans2 == 'T') dopCCD = .true. + if(ans3 == 'T') doDCD = .true. + if(ans4 == 'T') doCCSD = .true. + if(ans5 == 'T') doCCSDT = .true. ! Read weird CC methods @@ -88,11 +84,11 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & do_lCCD = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') do_drCCD = .true. - if(answer2 == 'T') do_rCCD = .true. - if(answer3 == 'T') do_crCCD = .true. - if(answer4 == 'T') do_lCCD = .true. + read(1,*) ans1,ans2,ans3,ans4 + if(ans1 == 'T') do_drCCD = .true. + if(ans2 == 'T') do_rCCD = .true. + if(ans3 == 'T') do_crCCD = .true. + if(ans4 == 'T') do_lCCD = .true. ! Read CI methods @@ -103,12 +99,12 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doFCI = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 - if(answer1 == 'T') doCIS = .true. - if(answer2 == 'T') doCIS_D = .true. - if(answer3 == 'T') doCID = .true. - if(answer4 == 'T') doCISD = .true. - if(answer5 == 'T') doFCI = .true. + read(1,*) ans1,ans2,ans3,ans4,ans5 + if(ans1 == 'T') doCIS = .true. + if(ans2 == 'T') doCIS_D = .true. + if(ans3 == 'T') doCID = .true. + if(ans4 == 'T') doCISD = .true. + if(ans5 == 'T') doFCI = .true. if(doCIS_D) doCIS = .true. ! Read RPA methods @@ -119,11 +115,11 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doppRPA = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') dophRPA = .true. - if(answer2 == 'T') dophRPAx = .true. - if(answer3 == 'T') docrRPA = .true. - if(answer4 == 'T') doppRPA = .true. + read(1,*) ans1,ans2,ans3,ans4 + if(ans1 == 'T') dophRPA = .true. + if(ans2 == 'T') dophRPAx = .true. + if(ans3 == 'T') docrRPA = .true. + if(ans4 == 'T') doppRPA = .true. ! Read Green's function methods @@ -134,12 +130,12 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doevGF3 = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 - if(answer1 == 'T') doG0F2 = .true. - if(answer2 == 'T') doevGF2 = .true. - if(answer3 == 'T') doqsGF2 = .true. - if(answer4 == 'T') doG0F3 = .true. - if(answer5 == 'T') doevGF3 = .true. + read(1,*) ans1,ans2,ans3,ans4,ans5 + if(ans1 == 'T') doG0F2 = .true. + if(ans2 == 'T') doevGF2 = .true. + if(ans3 == 'T') doqsGF2 = .true. + if(ans4 == 'T') doG0F3 = .true. + if(ans5 == 'T') doevGF3 = .true. ! Read GW methods @@ -151,13 +147,13 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doufGW = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 - if(answer1 == 'T') doG0W0 = .true. - if(answer2 == 'T') doevGW = .true. - if(answer3 == 'T') doqsGW = .true. - if(answer4 == 'T') doSRGqsGW = .true. - if(answer5 == 'T') doufG0W0 = .true. - if(answer6 == 'T') doufGW = .true. + read(1,*) ans1,ans2,ans3,ans4,ans5,ans6 + if(ans1 == 'T') doG0W0 = .true. + if(ans2 == 'T') doevGW = .true. + if(ans3 == 'T') doqsGW = .true. + if(ans4 == 'T') doSRGqsGW = .true. + if(ans5 == 'T') doufG0W0 = .true. + if(ans6 == 'T') doufGW = .true. ! Read GT methods @@ -169,13 +165,13 @@ subroutine read_methods(doRHF,doUHF,doROHF,doRMOM,doUMOM,doKS, & doqsGTeh = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 - if(answer1 == 'T') doG0T0pp = .true. - if(answer2 == 'T') doevGTpp = .true. - if(answer3 == 'T') doqsGTpp = .true. - if(answer4 == 'T') doG0T0eh = .true. - if(answer5 == 'T') doevGTeh = .true. - if(answer6 == 'T') doqsGTeh = .true. + read(1,*) ans1,ans2,ans3,ans4,ans5,ans6 + if(ans1 == 'T') doG0T0pp = .true. + if(ans2 == 'T') doevGTpp = .true. + if(ans3 == 'T') doqsGTpp = .true. + if(ans4 == 'T') doG0T0eh = .true. + if(ans5 == 'T') doevGTeh = .true. + if(ans6 == 'T') doqsGTeh = .true. ! Close file with geometry specification diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 6eed93f..246514b 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,11 +1,11 @@ -subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,max_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & - reg_MP, & - maxSCF_CC,thresh_CC,DIIS_CC,max_diis_CC, & - TDA,singlet,triplet,spin_conserved,spin_flip, & - maxSCF_GF,thresh_GF,DIIS_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_GT,thresh_GT,DIIS_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & - doACFDT,exchange_kernel,doXBS, & +subroutine read_options(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,dostab, & + reg_MP, & + maxSCF_CC,thresh_CC,max_diis_CC, & + TDA,spin_conserved,spin_flip, & + maxSCF_GF,thresh_GF,max_diis_GF,lin_GF,eta_GF,renorm_GF,reg_GF, & + maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, & + maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, & + doACFDT,exchange_kernel,doXBS, & dophBSE,dophBSE2,doppBSE,dBSE,dTDA) ! 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 double precision,intent(out) :: thresh_HF - logical,intent(out) :: DIIS_HF integer,intent(out) :: max_diis_HF integer,intent(out) :: guess_type - integer,intent(out) :: ortho_type - logical,intent(out) :: mix + double precision,intent(out) :: mix double precision,intent(out) :: level_shift 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 double precision,intent(out) :: thresh_CC - logical,intent(out) :: DIIS_CC integer,intent(out) :: max_diis_CC logical,intent(out) :: TDA - logical,intent(out) :: singlet - logical,intent(out) :: triplet logical,intent(out) :: spin_conserved logical,intent(out) :: spin_flip integer,intent(out) :: maxSCF_GF double precision,intent(out) :: thresh_GF - logical,intent(out) :: DIIS_GF integer,intent(out) :: max_diis_GF logical,intent(out) :: lin_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 double precision,intent(out) :: thresh_GW - logical,intent(out) :: DIIS_GW integer,intent(out) :: max_diis_GW logical,intent(out) :: TDA_W 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 double precision,intent(out) :: thresh_GT - logical,intent(out) :: DIIS_GT integer,intent(out) :: max_diis_GT logical,intent(out) :: TDA_T 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 thresh_HF = 1d-6 - DIIS_HF = .false. - max_diis_HF = 5 + max_diis_HF = 1 guess_type = 1 - ortho_type = 1 - mix = .false. + mix = 0d0 level_shift = 0d0 dostab = .false. 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(ans2 == 'T') mix = .true. - if(ans3 == 'T') dostab = .true. - - if(.not.DIIS_HF) max_diis_HF = 1 + if(ans1 == 'T') dostab = .true. ! 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 thresh_CC = 1d-5 - DIIS_CC = .false. - max_diis_CC = 5 + max_diis_CC = 1 read(1,*) - read(1,*) maxSCF_CC,thresh_CC,ans1,max_diis_CC - - if(ans1 == 'T') DIIS_CC = .true. - - if(.not.DIIS_CC) max_diis_CC = 1 + read(1,*) maxSCF_CC,thresh_CC,max_diis_CC ! Read excited state options TDA = .false. - singlet = .false. - triplet = .false. spin_conserved = .false. spin_flip = .false. read(1,*) - read(1,*) ans1,ans2,ans3,ans4,ans5 + read(1,*) ans1,ans2,ans3 if(ans1 == 'T') TDA = .true. - if(ans2 == 'T') singlet = .true. - if(ans3 == 'T') triplet = .true. - if(ans4 == 'T') spin_conserved = .true. - if(ans5 == 'T') spin_flip = .true. + if(ans2 == 'T') spin_conserved = .true. + if(ans3 == 'T') spin_flip = .true. ! Read GF options maxSCF_GF = 64 thresh_GF = 1d-5 - DIIS_GF = .false. - max_diis_GF = 5 + max_diis_GF = 1 lin_GF = .false. eta_GF = 0d0 renorm_GF = 0 reg_GF = .false. 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(ans2 == 'T') lin_GF = .true. - if(ans3 == 'T') reg_GF = .true. - if(.not.DIIS_GF) max_diis_GF = 1 + if(ans1 == 'T') lin_GF = .true. + if(ans2 == 'T') reg_GF = .true. ! Read GW options maxSCF_GW = 64 thresh_GW = 1d-5 - DIIS_GW = .false. - max_diis_GW = 5 + max_diis_GW = 1 lin_GW = .false. eta_GW = 0d0 reg_GW = .false. TDA_W = .false. 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(ans2 == 'T') lin_GW = .true. - if(ans3 == 'T') TDA_W = .true. - if(ans4 == 'T') reg_GW = .true. - if(.not.DIIS_GW) max_diis_GW = 1 + if(ans1 == 'T') lin_GW = .true. + if(ans2 == 'T') TDA_W = .true. + if(ans3 == 'T') reg_GW = .true. ! Read GT options maxSCF_GT = 64 thresh_GT = 1d-5 - DIIS_GT = .false. - max_diis_GT = 5 + max_diis_GT = 1 lin_GT = .false. eta_GT = 0d0 reg_GT = .false. TDA_T = .false. 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(ans2 == 'T') lin_GT = .true. - if(ans3 == 'T') TDA_T = .true. - if(ans4 == 'T') reg_GT = .true. - if(.not.DIIS_GT) max_diis_GT = 1 + if(ans1 == 'T') lin_GT = .true. + if(ans2 == 'T') TDA_T = .true. + if(ans3 == 'T') reg_GT = .true. ! Options for adiabatic connection diff --git a/src/RPA/RRPA.f90 b/src/RPA/RRPA.f90 new file mode 100644 index 0000000..4c36d77 --- /dev/null +++ b/src/RPA/RRPA.f90 @@ -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 diff --git a/src/RPA/RPA.f90 b/src/RPA/URPA.f90 similarity index 62% rename from src/RPA/RPA.f90 rename to src/RPA/URPA.f90 index a67550d..da77aa8 100644 --- a/src/RPA/RPA.f90 +++ b/src/RPA/URPA.f90 @@ -1,7 +1,5 @@ -subroutine 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,ERI_aaaa,ERI_aabb,ERI_bbbb, & - dipole_int,dipole_int_aa,dipole_int_bb,epsHF,cHF,S) +subroutine 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) ! Random-phase approximation module @@ -14,13 +12,10 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, logical :: dophRPAx logical :: docrRPA logical :: doppRPA - logical :: unrestricted logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel - logical,intent(in) :: singlet - logical,intent(in) :: triplet logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip 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) :: cHF(nBas,nBas,nspin) 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_aabb(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_bb(nBas,nBas,ncart) @@ -53,12 +46,8 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(dophRPA) then 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, & - 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 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) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -74,12 +63,8 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(dophRPAx) then 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, & - 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 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) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -95,11 +80,7 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(docrRPA) then call wall_time(start_RPA) - if(unrestricted) then - 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 + write(*,*) 'Unrestricted version of crRPA not yet implemented! Sorry.' call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -115,11 +96,7 @@ subroutine RPA(dophRPA,dophRPAx,docrRPA,doppRPA,unrestricted, if(doppRPA) then 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) - else - call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,EHF,ERI,dipole_int,epsHF) - end if + call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EHF,ERI_aaaa,ERI_aabb,ERI_bbbb,epsHF) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA diff --git a/src/utils/orthogonalization_matrix.f90 b/src/utils/orthogonalization_matrix.f90 index 5c78601..42ecb67 100644 --- a/src/utils/orthogonalization_matrix.f90 +++ b/src/utils/orthogonalization_matrix.f90 @@ -1,4 +1,4 @@ -subroutine orthogonalization_matrix(ortho_type,nBas,S,X) +subroutine orthogonalization_matrix(nBas,S,X) ! Compute the orthogonalization matrix X @@ -6,7 +6,7 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X) ! Input variables - integer,intent(in) :: nBas,ortho_type + integer,intent(in) :: nBas double precision,intent(in) :: S(nBas,nBas) ! Local variables @@ -14,6 +14,7 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X) logical :: debug double precision,allocatable :: UVec(:,:),Uval(:) double precision,parameter :: thresh = 1d-6 + integer,parameter :: ortho_type = 1 integer :: i diff --git a/src/utils/read_molecule.f90 b/src/utils/read_molecule.f90 index 3d29227..4c018f3 100644 --- a/src/utils/read_molecule.f90 +++ b/src/utils/read_molecule.f90 @@ -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 @@ -6,19 +6,18 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR) 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 integer :: nCore 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(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(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 @@ -35,9 +34,8 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR) end if - nO(:) = nEl(:) nC(:) = nCore/2 - nR(:) = nRyd + nR(:) = nRyd/2 ! Print results @@ -46,9 +44,9 @@ subroutine read_molecule(nNuc,nEl,nO,nC,nR) write(*,'(A28)') '----------------------' write(*,*) write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nEl(1) - write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nEl(2) - write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nEl(:)) + write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nO(1) + write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nO(2) + write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nO(:)) write(*,'(A28)') '----------------------' write(*,*) write(*,'(A28)') '----------------------'