From 28d85ed2043ae9ad09b14774901f96d2b5e7dd15 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 28 Sep 2020 22:24:02 +0200 Subject: [PATCH] unrestricted f --- input/basis | 44 +++--- input/methods | 6 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/options | 4 +- src/QuAcK/URPAx.f90 | 2 +- src/QuAcK/print_UHF.f90 | 1 + .../print_unrestricted_transition_vectors.f90 | 132 ++++++++++++++++++ 8 files changed, 171 insertions(+), 28 deletions(-) create mode 100644 src/QuAcK/print_unrestricted_transition_vectors.f90 diff --git a/input/basis b/input/basis index fb05e68..b2b2293 100644 --- a/input/basis +++ b/input/basis @@ -1,18 +1,30 @@ -1 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 +1 6 +S 8 + 1 1469.0000000 0.0007660 + 2 220.5000000 0.0058920 + 3 50.2600000 0.0296710 + 4 14.2400000 0.1091800 + 5 4.5810000 0.2827890 + 6 1.5800000 0.4531230 + 7 0.5640000 0.2747740 + 8 0.0734500 0.0097510 +S 8 + 1 1469.0000000 -0.0001200 + 2 220.5000000 -0.0009230 + 3 50.2600000 -0.0046890 + 4 14.2400000 -0.0176820 + 5 4.5810000 -0.0489020 + 6 1.5800000 -0.0960090 + 7 0.5640000 -0.1363800 + 8 0.0734500 0.5751020 S 1 - 1 0.1220000 1.0000000 + 1 0.0280500 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 P 1 - 1 0.7270000 1.0000000 -2 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 -S 1 - 1 0.1220000 1.0000000 -P 1 - 1 0.7270000 1.0000000 + 1 0.0240300 1.0000000 +D 1 + 1 0.1239000 1.0000000 + diff --git a/input/methods b/input/methods index 0d456f5..adb472a 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - T F F + F T F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) @@ -9,11 +9,11 @@ # CIS* CID CISD F F F # RPA* RPAx* ppRPA - F F F + F T 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/molecule b/input/molecule index 7225285..058d6dd 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 1 2 1 0 0 # Znuc x y z - H 0. 0. -0.7 - H 0. 0. 0.7 + Li 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index ac9420c..c9a5a65 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - H 0.0000000000 0.0000000000 -0.3704240743 - H 0.0000000000 0.0000000000 0.3704240743 + Li 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index 1d00d01..1a82a15 100644 --- a/input/options +++ b/input/options @@ -1,11 +1,11 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type - 64 0.00001 T 5 1 1 + 64 0.0000001 T 5 1 1 # MP: # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T T F + T T T F 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/URPAx.f90 b/src/QuAcK/URPAx.f90 index 9a3c510..ad58315 100644 --- a/src/QuAcK/URPAx.f90 +++ b/src/QuAcK/URPAx.f90 @@ -78,7 +78,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n 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(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_sc,dipole_int,Omega_sc,XpY_sc,XmY_sc) deallocate(Omega_sc,XpY_sc,XmY_sc) diff --git a/src/QuAcK/print_UHF.f90 b/src/QuAcK/print_UHF.f90 index 0ad8c89..640d023 100644 --- a/src/QuAcK/print_UHF.f90 +++ b/src/QuAcK/print_UHF.f90 @@ -84,6 +84,7 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) write(*,'(A50)') 'UHF spin-up orbital coefficients ' write(*,'(A50)') '-----------------------------------------' call matout(nBas,nBas,c(:,:,1)) + write(*,*) write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') 'UHF spin-down orbital coefficients ' write(*,'(A50)') '-----------------------------------------' diff --git a/src/QuAcK/print_unrestricted_transition_vectors.f90 b/src/QuAcK/print_unrestricted_transition_vectors.f90 new file mode 100644 index 0000000..2d6e8bf --- /dev/null +++ b/src/QuAcK/print_unrestricted_transition_vectors.f90 @@ -0,0 +1,132 @@ +subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSt,dipole_int,Omega,XpY,XmY) + +! Print transition vectors for linear response calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: spin_allowed + 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) + integer,intent(in) :: nSt + double precision :: dipole_int(nBas,nBas,ncart,nspin) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: XpY(nSt,nSt) + double precision,intent(in) :: XmY(nSt,nSt) + +! Local variables + + integer :: ia,jb,i,j,a,b + integer :: ixyz + integer :: ispin + integer,parameter :: maxS = 10 + double precision :: norm + double precision,parameter :: thres_vec = 0.1d0 + double precision,allocatable :: X(:) + double precision,allocatable :: Y(:) + double precision,allocatable :: f(:,:) + double precision,allocatable :: os(:) + +! Memory allocation + + allocate(X(nSt),Y(nSt),f(nSt,ncart),os(nSt)) + +! Compute dipole moments and oscillator strengths + + + f(:,:) = 0d0 + if(spin_allowed) then + + do ispin=1,nspin + do ia=1,nSt + do ixyz=1,ncart + jb = 0 + do j=nC(ispin)+1,nO(ispin) + do b=nO(ispin)+1,nBas-nR(ispin) + jb = jb + 1 + f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz,ispin)*XpY(ia,jb) + end do + end do + end do + end do + end do + + write(*,*) '----------------' + write(*,*) ' Dipole moments ' + write(*,*) '----------------' + call matout(nSt,ncart,f(:,:)) + write(*,*) + + do ia=1,nSt + os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2) + end do + + write(*,*) '----------------------' + write(*,*) ' Oscillator strengths ' + write(*,*) '----------------------' + call matout(nSt,1,os(:)) + write(*,*) + + end if + +! Print details about excitations + + do ia=1,min(nSt,maxS) + + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + print*,'---------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')' + print*,'---------------------------------------------' + + ! Spin-up transitions + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0) + end do + end do + + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0) + end do + end do + write(*,*) + + ! Spin-down transitions + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0) + end do + end do + + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0) + end do + end do + write(*,*) + + end do + + write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) + write(*,*) + +end subroutine print_unrestricted_transition_vectors