From a611ee7442ed0b8edba66f40022ffa9d7ba32a95 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 24 Sep 2020 16:39:15 +0200 Subject: [PATCH] spin flip --- input/methods | 4 ++-- input/options | 2 +- src/QuAcK/QuAcK.f90 | 18 +++++++++++++----- src/QuAcK/RPAx.f90 | 9 +++++---- src/QuAcK/UCIS.f90 | 4 ++-- src/QuAcK/URPAx.f90 | 9 +++++---- src/QuAcK/UdRPA.f90 | 9 +++++---- src/QuAcK/dRPA.f90 | 10 +++++----- .../unrestricted_linear_response_A_matrix.f90 | 2 +- 9 files changed, 39 insertions(+), 28 deletions(-) diff --git a/input/methods b/input/methods index ba17592..70cc161 100644 --- a/input/methods +++ b/input/methods @@ -7,13 +7,13 @@ # drCCD rCCD lCCD pCCD F F F F # CIS CID CISD - F F F + T F F # RPA RPAx ppRPA F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 55538aa..cb56d2c 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T T T + T T T T F # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 3bc4577..4a63fe3 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -586,7 +586,15 @@ program QuAcK if(doCIS) then call cpu_time(start_CIS) - call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF) + if(unrestricted) then + + call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,eHF) + + else + + call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF) + + end if call cpu_time(end_CIS) t_CIS = end_CIS - start_CIS @@ -636,12 +644,12 @@ program QuAcK call cpu_time(start_RPA) if(unrestricted) then - call UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + call UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,eHF) else - call dRPA(doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) end if call cpu_time(end_RPA) @@ -661,12 +669,12 @@ program QuAcK call cpu_time(start_RPAx) if(unrestricted) then - call URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & + call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,eHF) else - call RPAx(doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,0d0,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) end if call cpu_time(end_RPAx) diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 index ca9f709..37f3f25 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -1,4 +1,4 @@ -subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform random phase approximation calculation with exchange (aka TDHF) @@ -8,6 +8,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ! Input variables + logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet @@ -58,7 +59,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ispin = 1 - call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -71,7 +72,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, ispin = 2 - call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -103,7 +104,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, write(*,*) '-------------------------------------------------------' write(*,*) - call ACFDT(exchange_kernel,.false.,.false.,.false.,.false.,.false.,singlet,triplet,eta, & + call ACFDT(exchange_kernel,.false.,.false.,.false.,TDA,.false.,singlet,triplet,eta, & nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) write(*,*) diff --git a/src/QuAcK/UCIS.f90 b/src/QuAcK/UCIS.f90 index d27a61b..9db0626 100644 --- a/src/QuAcK/UCIS.f90 +++ b/src/QuAcK/UCIS.f90 @@ -64,7 +64,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E allocate(A_sc(nS_sc,nS_sc),Omega_sc(nS_sc)) - call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,lambda,eHF, & + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,lambda,eHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sc) if(dump_matrix) then @@ -102,7 +102,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf)) - call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,lambda,eHF, & + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,lambda,eHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf) if(dump_matrix) then diff --git a/src/QuAcK/URPAx.f90 b/src/QuAcK/URPAx.f90 index f351ca6..ff13fad 100644 --- a/src/QuAcK/URPAx.f90 +++ b/src/QuAcK/URPAx.f90 @@ -1,4 +1,4 @@ -subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & +subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,e) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -9,7 +9,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO ! Input variables - double precision,intent(in) :: eta + logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: spin_conserved @@ -20,6 +20,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) + double precision,intent(in) :: eta double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: e(nBas,nspin) @@ -73,7 +74,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPAx ',5,nS_sc,Omega_sc) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -96,7 +97,7 @@ subroutine URPAx(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & + call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPAx ',6,nS_sf,Omega_sf) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/QuAcK/UdRPA.f90 b/src/QuAcK/UdRPA.f90 index 9fff2be..0169441 100644 --- a/src/QuAcK/UdRPA.f90 +++ b/src/QuAcK/UdRPA.f90 @@ -1,4 +1,4 @@ -subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & +subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,e) ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism @@ -9,7 +9,7 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO ! Input variables - double precision,intent(in) :: eta + logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: spin_conserved @@ -20,6 +20,7 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO integer,intent(in) :: nV(nspin) integer,intent(in) :: nR(nspin) integer,intent(in) :: nS(nspin) + double precision,intent(in) :: eta double precision,intent(in) :: ENuc double precision,intent(in) :: EUHF double precision,intent(in) :: e(nBas,nspin) @@ -73,7 +74,7 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sc(nS_sc),XpY_sc(nS_sc,nS_sc),XmY_sc(nS_sc,nS_sc)) - call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & + call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc) call print_excitation('URPA ',5,nS_sc,Omega_sc) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -95,7 +96,7 @@ subroutine UdRPA(doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & + call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf) call print_excitation('URPA ',6,nS_sf,Omega_sf) ! call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/QuAcK/dRPA.f90 b/src/QuAcK/dRPA.f90 index d17ed51..3985af2 100644 --- a/src/QuAcK/dRPA.f90 +++ b/src/QuAcK/dRPA.f90 @@ -1,5 +1,4 @@ -subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine dRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform a direct random phase approximation calculation @@ -9,6 +8,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, & ! Input variables + logical,intent(in) :: TDA logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: singlet @@ -59,7 +59,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, & ispin = 1 - call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -72,7 +72,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, & ispin = 2 - call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -105,7 +105,7 @@ subroutine dRPA(doACFDT,exchange_kernel,singlet,triplet,eta, & write(*,*) '------------------------------------------------------' write(*,*) - call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,.false.,singlet,triplet,eta, & + call ACFDT(exchange_kernel,.false.,.true.,.false.,TDA,.false.,singlet,triplet,eta, & nBas,nC,nO,nV,nR,nS,ERI,eHF,eHF,EcAC) if(exchange_kernel) then diff --git a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 index 6c282e5..2f2d1b8 100644 --- a/src/QuAcK/unrestricted_linear_response_A_matrix.f90 +++ b/src/QuAcK/unrestricted_linear_response_A_matrix.f90 @@ -162,7 +162,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa jb = jb + 1 A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & - - (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j) + - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a) end do end do