From 769f4e922116ac1eeb803a16a07ad15e14acecc2 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 5 Jun 2020 22:34:32 +0200 Subject: [PATCH] (1) --- input/basis | 127 ++++++++++++------------- input/methods | 4 +- src/QuAcK/QuAcK.f90 | 38 ++++---- src/QuAcK/RPAx.f90 | 3 + src/QuAcK/print_transition_vectors.f90 | 32 ++++--- 5 files changed, 104 insertions(+), 100 deletions(-) diff --git a/input/basis b/input/basis index bbe0bfe..80356bb 100644 --- a/input/basis +++ b/input/basis @@ -1,71 +1,64 @@ -1 9 -S 8 - 1 9046.0000000 0.0007000 - 2 1357.0000000 0.0053890 - 3 309.3000000 0.0274060 - 4 87.7300000 0.1032070 - 5 28.5600000 0.2787230 - 6 10.2100000 0.4485400 - 7 3.8380000 0.2782380 - 8 0.7466000 0.0154400 -S 8 - 1 9046.0000000 -0.0001530 - 2 1357.0000000 -0.0012080 - 3 309.3000000 -0.0059920 - 4 87.7300000 -0.0245440 - 5 28.5600000 -0.0674590 - 6 10.2100000 -0.1580780 - 7 3.8380000 -0.1218310 - 8 0.7466000 0.5490030 +1 6 +S 9 +1 9.046000E+03 7.000000E-04 +2 1.357000E+03 5.389000E-03 +3 3.093000E+02 2.740600E-02 +4 8.773000E+01 1.032070E-01 +5 2.856000E+01 2.787230E-01 +6 1.021000E+01 4.485400E-01 +7 3.838000E+00 2.782380E-01 +8 7.466000E-01 1.544000E-02 +9 2.248000E-01 -2.864000E-03 +S 9 +1 9.046000E+03 -1.530000E-04 +2 1.357000E+03 -1.208000E-03 +3 3.093000E+02 -5.992000E-03 +4 8.773000E+01 -2.454400E-02 +5 2.856000E+01 -6.745900E-02 +6 1.021000E+01 -1.580780E-01 +7 3.838000E+00 -1.218310E-01 +8 7.466000E-01 5.490030E-01 +9 2.248000E-01 5.788150E-01 S 1 - 1 0.2248000 1.0000000 +1 2.248000E-01 1.000000E+00 +P 4 +1 1.355000E+01 3.991900E-02 +2 2.917000E+00 2.171690E-01 +3 7.973000E-01 5.103190E-01 +4 2.185000E-01 4.622140E-01 +P 1 +1 2.185000E-01 1.000000E+00 +D 1 +1 8.170000E-01 1.0000000 +2 6 +S 9 +1 9.046000E+03 7.000000E-04 +2 1.357000E+03 5.389000E-03 +3 3.093000E+02 2.740600E-02 +4 8.773000E+01 1.032070E-01 +5 2.856000E+01 2.787230E-01 +6 1.021000E+01 4.485400E-01 +7 3.838000E+00 2.782380E-01 +8 7.466000E-01 1.544000E-02 +9 2.248000E-01 -2.864000E-03 +S 9 +1 9.046000E+03 -1.530000E-04 +2 1.357000E+03 -1.208000E-03 +3 3.093000E+02 -5.992000E-03 +4 8.773000E+01 -2.454400E-02 +5 2.856000E+01 -6.745900E-02 +6 1.021000E+01 -1.580780E-01 +7 3.838000E+00 -1.218310E-01 +8 7.466000E-01 5.490030E-01 +9 2.248000E-01 5.788150E-01 S 1 - 1 0.0612400 1.0000000 -P 3 - 1 13.5500000 0.0399190 - 2 2.9170000 0.2171690 - 3 0.7973000 0.5103190 +1 2.248000E-01 1.000000E+00 +P 4 +1 1.355000E+01 3.991900E-02 +2 2.917000E+00 2.171690E-01 +3 7.973000E-01 5.103190E-01 +4 2.185000E-01 4.622140E-01 P 1 - 1 0.2185000 1.0000000 -P 1 - 1 0.0561100 1.0000000 +1 2.185000E-01 1.000000E+00 D 1 - 1 0.8170000 1.0000000 -D 1 - 1 0.2300000 1.0000000 -2 9 -S 8 - 1 9046.0000000 0.0007000 - 2 1357.0000000 0.0053890 - 3 309.3000000 0.0274060 - 4 87.7300000 0.1032070 - 5 28.5600000 0.2787230 - 6 10.2100000 0.4485400 - 7 3.8380000 0.2782380 - 8 0.7466000 0.0154400 -S 8 - 1 9046.0000000 -0.0001530 - 2 1357.0000000 -0.0012080 - 3 309.3000000 -0.0059920 - 4 87.7300000 -0.0245440 - 5 28.5600000 -0.0674590 - 6 10.2100000 -0.1580780 - 7 3.8380000 -0.1218310 - 8 0.7466000 0.5490030 -S 1 - 1 0.2248000 1.0000000 -S 1 - 1 0.0612400 1.0000000 -P 3 - 1 13.5500000 0.0399190 - 2 2.9170000 0.2171690 - 3 0.7973000 0.5103190 -P 1 - 1 0.2185000 1.0000000 -P 1 - 1 0.0561100 1.0000000 -D 1 - 1 0.8170000 1.0000000 -D 1 - 1 0.2300000 1.0000000 - +1 8.170000E-01 1.0000000 diff --git a/input/methods b/input/methods index 63c03d4..aba2ee4 100644 --- a/input/methods +++ b/input/methods @@ -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/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 718f0b6..023058a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -383,7 +383,7 @@ program QuAcK if(doCCD) then call cpu_time(start_CCD) - call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & + call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) @@ -402,7 +402,7 @@ program QuAcK if(doCCSD) then call cpu_time(start_CCSD) - call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC(1),nO(1),nV(1),nR(1), & + call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCSD) @@ -419,7 +419,7 @@ program QuAcK if(do_drCCD) then call cpu_time(start_CCD) - call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & + call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) @@ -436,7 +436,7 @@ program QuAcK if(do_rCCD) then call cpu_time(start_CCD) - call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & + call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) @@ -453,7 +453,7 @@ program QuAcK if(do_lCCD) then call cpu_time(start_CCD) - call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & + call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) @@ -470,7 +470,7 @@ program QuAcK if(do_pCCD) then call cpu_time(start_CCD) - call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & + call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) @@ -570,7 +570,7 @@ program QuAcK call cpu_time(start_ppRPA) call ppRPA(singlet_manifold,triplet_manifold, & - nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO,eHF) + nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_ppRPA) t_ppRPA = end_ppRPA - start_ppRPA @@ -587,7 +587,7 @@ program QuAcK ! call cpu_time(start_ADC) ! call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, & -! nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO) +! nBas,nC,nO,nV,nR,eHF,ERI_MO) ! call cpu_time(end_ADC) ! t_ADC = end_ADC - start_ADC @@ -604,7 +604,7 @@ program QuAcK call cpu_time(start_GF2) call G0F2(BSE_GF,TDA_GF,singlet_manifold,triplet_manifold,linGF, & - eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) + eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -621,7 +621,7 @@ program QuAcK call cpu_time(start_GF2) call evGF2(BSE_GF,TDA_GF,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, & - eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) + eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -637,7 +637,7 @@ program QuAcK if(doG0F3) then call cpu_time(start_GF3) - call G0F3(renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) + call G0F3(renormGF,nBas,nC,nO,nV,nR,ERI_MO,eHF) call cpu_time(end_GF3) t_GF3 = end_GF3 - start_GF3 @@ -653,7 +653,7 @@ program QuAcK if(doevGF3) then call cpu_time(start_GF3) - call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) + call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC,nO,nV,nR,ERI_MO,eHF) call cpu_time(end_GF3) t_GF3 = end_GF3 - start_GF3 @@ -673,7 +673,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0 @@ -691,7 +691,7 @@ program QuAcK call cpu_time(start_evGW) call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) t_evGW = end_evGW - start_evGW @@ -709,7 +709,7 @@ program QuAcK call cpu_time(start_qsGW) call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) call cpu_time(end_qsGW) t_qsGW = end_qsGW - start_qsGW @@ -729,7 +729,7 @@ program QuAcK call cpu_time(start_G0T0) call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_G0T0) t_G0T0 = end_G0T0 - start_G0T0 @@ -747,7 +747,7 @@ program QuAcK call cpu_time(start_evGT) call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & BSE_GW,TDA_GW,singlet_manifold,triplet_manifold,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_evGT) t_evGT = end_evGT - start_evGT @@ -864,7 +864,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0 @@ -878,7 +878,7 @@ program QuAcK call cpu_time(start_G0T0) call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, & singlet_manifold,triplet_manifold,linGW,eta_GW, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) call cpu_time(end_G0T0) t_G0T0 = end_G0T0 - start_G0T0 diff --git a/src/QuAcK/RPAx.f90 b/src/QuAcK/RPAx.f90 index c57d638..10d2f58 100644 --- a/src/QuAcK/RPAx.f90 +++ b/src/QuAcK/RPAx.f90 @@ -62,6 +62,8 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + endif @@ -74,6 +76,7 @@ subroutine RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & call linear_response(ispin,.false.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx ',ispin,nS,Omega(:,ispin)) + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif diff --git a/src/QuAcK/print_transition_vectors.f90 b/src/QuAcK/print_transition_vectors.f90 index 0dfaf8c..ab3e093 100644 --- a/src/QuAcK/print_transition_vectors.f90 +++ b/src/QuAcK/print_transition_vectors.f90 @@ -22,6 +22,7 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) integer :: ia,jb,i,j,a,b integer,parameter :: maxS = 10 double precision,parameter :: thres_vec = 0.1d0 + double precision :: norm double precision,allocatable :: X(:) double precision,allocatable :: Y(:) @@ -32,31 +33,38 @@ subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) write(*,*) do ia=1,min(nS,maxS) - X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) - Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + norm = 0d0 + do jb=1,nS + norm = norm + X(jb)*X(jb) + end do + norm = sqrt(norm) - print*,'--------------------------------' - write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV' - print*,'--------------------------------' - print*,'Resonant vectors' - print*,'--------------------------------' + print*,'--------------------------------' + write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV' + print*,'--------------------------------' jb = 0 do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',X(jb) + if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0) end do end do - print*,'--------------------------------' - print*,'Anti-resonant vectors' - print*,'--------------------------------' + norm = 0d0 + do jb=1,nS + norm = norm + Y(jb)*Y(jb) + end do + norm = sqrt(norm) + jb = 0 do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',Y(jb) + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0) end do end do write(*,*)