From ac66ae8ea214e8bcb1e8fefba708c556bde1fc65 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 30 Sep 2020 15:02:48 +0200 Subject: [PATCH] corrected bug in spin flip --- examples/molecule.Be | 2 +- input/methods | 4 +- src/CI/UCIS.f90 | 12 ++--- src/MBPT/unrestricted_Bethe_Salpeter.f90 | 2 +- src/QuAcK/QuAcK.f90 | 4 +- src/RPA/URPAx.f90 | 2 +- src/RPA/UdRPA.f90 | 6 +-- src/RPA/print_excitation.f90 | 45 ------------------- .../unrestricted_linear_response_A_matrix.f90 | 3 +- 9 files changed, 19 insertions(+), 61 deletions(-) delete mode 100644 src/RPA/print_excitation.f90 diff --git a/examples/molecule.Be b/examples/molecule.Be index 6a6f6d1..7f85101 100644 --- a/examples/molecule.Be +++ b/examples/molecule.Be @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 2 2 0 0 + 1 3 1 0 0 # Znuc x y z Be 0.0 0.0 0.0 diff --git a/input/methods b/input/methods index adb472a..bb07f30 100644 --- a/input/methods +++ b/input/methods @@ -7,9 +7,9 @@ # drCCD rCCD lCCD pCCD F F F F # CIS* CID CISD - F F F + T F F # RPA* RPAx* ppRPA - F T F + F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0* evGW* qsGW diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index e75fc5a..5a8e7bd 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -1,4 +1,5 @@ -subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,eHF) +subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & + dipole_int_aa,dipole_int_bb,eHF) ! Perform configuration interaction single calculation` @@ -20,7 +21,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart,nspin) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart,nspin) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart,nspin) ! Local variables @@ -102,10 +104,10 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E nS_sf = nS_ab + nS_ba 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,lambda,eHF, & + + call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf) - + if(dump_matrix) then print*,'CIS matrix (spin-conserved transitions)' call matout(nS_sf,nS_sf,A_sf) diff --git a/src/MBPT/unrestricted_Bethe_Salpeter.f90 b/src/MBPT/unrestricted_Bethe_Salpeter.f90 index 0ce2590..7e4dc02 100644 --- a/src/MBPT/unrestricted_Bethe_Salpeter.f90 +++ b/src/MBPT/unrestricted_Bethe_Salpeter.f90 @@ -142,7 +142,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, !------------------------------------------------- if(dBSE) & - call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sf,nS_sc, & + call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, & eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb, & OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 4a906c9..f16eab9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -617,7 +617,7 @@ program QuAcK 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,dipole_int,eHF) + ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) else @@ -674,7 +674,7 @@ program QuAcK if(unrestricted) then 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,dipole_int,eHF) + ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF) else diff --git a/src/RPA/URPAx.f90 b/src/RPA/URPAx.f90 index 12172c2..2f21783 100644 --- a/src/RPA/URPAx.f90 +++ b/src/RPA/URPAx.f90 @@ -108,7 +108,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.false.,TDA,.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_ab,nS_ba,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_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & diff --git a/src/RPA/UdRPA.f90 b/src/RPA/UdRPA.f90 index 4efddc4..99d94de 100644 --- a/src/RPA/UdRPA.f90 +++ b/src/RPA/UdRPA.f90 @@ -106,11 +106,11 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf)) - call unrestricted_linear_response(ispin,.true.,TDA,.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_ab,nS_ba,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)) - + call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & + Omega_sf,XpY_sf,XmY_sf) endif diff --git a/src/RPA/print_excitation.f90 b/src/RPA/print_excitation.f90 deleted file mode 100644 index c2bbdb4..0000000 --- a/src/RPA/print_excitation.f90 +++ /dev/null @@ -1,45 +0,0 @@ -subroutine print_excitation(method,ispin,nS,Omega) - -! Print excitation energies for a given spin manifold - - implicit none - include 'parameters.h' - -! Input variables - - character*12,intent(in) :: method - integer,intent(in) :: ispin,nS - double precision,intent(in) :: Omega(nS) - -! Local variables - - character*14 :: spin_manifold - integer,parameter :: maxS = 32 - integer :: ia - - if(ispin == 1) spin_manifold = 'singlet' - if(ispin == 2) spin_manifold = 'triplet' - if(ispin == 3) spin_manifold = 'alpha-beta' - if(ispin == 4) spin_manifold = 'alpha-alpha' - if(ispin == 5) spin_manifold = 'spin-conserved' - if(ispin == 6) spin_manifold = 'spin-flip' - - write(*,*) - write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A14,A14,A14,A9)') method,' calculation: ',spin_manifold,' manifold' - write(*,*)'-------------------------------------------------------------' - write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') & - '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' - write(*,*)'-------------------------------------------------------------' - - do ia=1,min(nS,maxS) - write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & - '|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|' - enddo - - write(*,*)'-------------------------------------------------------------' - write(*,*) - -end subroutine print_excitation - - diff --git a/src/RPA/unrestricted_linear_response_A_matrix.f90 b/src/RPA/unrestricted_linear_response_A_matrix.f90 index 2f2d1b8..8754667 100644 --- a/src/RPA/unrestricted_linear_response_A_matrix.f90 +++ b/src/RPA/unrestricted_linear_response_A_matrix.f90 @@ -162,7 +162,8 @@ 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,j,i,a) + - (1d0 - delta_dRPA)*lambda*ERI_abab(j,a,i,b) +! print*,'ia,jb,A(ia,jb) = ',i,a,j,b,ia,jb,A_lr(ia,jb) end do end do