diff --git a/input/basis b/input/basis index 4b80620..894420f 100644 --- a/input/basis +++ b/input/basis @@ -1,46 +1,98 @@ -1 10 -S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 +1 15 +S 9 + 1 33980.0000000 0.0000910 + 2 5089.0000000 0.0007040 + 3 1157.0000000 0.0036930 + 4 326.6000000 0.0153600 + 5 106.1000000 0.0529290 + 6 38.1100000 0.1470430 + 7 14.7500000 0.3056310 + 8 6.0350000 0.3993450 + 9 2.5300000 0.2170510 +S 9 + 1 33980.0000000 -0.0000190 + 2 5089.0000000 -0.0001510 + 3 1157.0000000 -0.0007850 + 4 326.6000000 -0.0033240 + 5 106.1000000 -0.0115120 + 6 38.1100000 -0.0341600 + 7 14.7500000 -0.0771730 + 8 6.0350000 -0.1414930 + 9 2.5300000 -0.1180190 S 1 - 1 0.7977000 1.0000000 + 1 0.7355000 1.0000000 S 1 - 1 0.2581000 1.0000000 + 1 0.2905000 1.0000000 S 1 - 1 0.0898900 1.0000000 + 1 0.1111000 1.0000000 +P 3 + 1 34.5100000 0.0053780 + 2 7.9150000 0.0361320 + 3 2.3680000 0.1424930 P 1 - 1 2.2920000 1.0000000 + 1 0.8132000 1.0000000 P 1 - 1 0.8380000 1.0000000 + 1 0.2890000 1.0000000 P 1 - 1 0.2920000 1.0000000 + 1 0.1007000 1.0000000 D 1 - 1 2.0620000 1.0000000 + 1 1.8480000 1.0000000 D 1 - 1 0.6620000 1.0000000 + 1 0.6490000 1.0000000 +D 1 + 1 0.2280000 1.0000000 F 1 - 1 1.3970000 1.0000000 -2 10 -S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 -S 1 - 1 0.7977000 1.0000000 -S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 -P 1 - 1 2.2920000 1.0000000 -P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 -D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 + 1 1.4190000 1.0000000 F 1 - 1 1.3970000 1.0000000 + 1 0.4850000 1.0000000 +G 1 + 1 1.0110000 1.0000000 +2 15 +S 9 + 1 61420.0000000 0.0000900 + 2 9199.0000000 0.0006980 + 3 2091.0000000 0.0036640 + 4 590.9000000 0.0152180 + 5 192.3000000 0.0524230 + 6 69.3200000 0.1459210 + 7 26.9700000 0.3052580 + 8 11.1000000 0.3985080 + 9 4.6820000 0.2169800 +S 9 + 1 61420.0000000 -0.0000200 + 2 9199.0000000 -0.0001590 + 3 2091.0000000 -0.0008290 + 4 590.9000000 -0.0035080 + 5 192.3000000 -0.0121560 + 6 69.3200000 -0.0362610 + 7 26.9700000 -0.0829920 + 8 11.1000000 -0.1520900 + 9 4.6820000 -0.1153310 +S 1 + 1 1.4280000 1.0000000 +S 1 + 1 0.5547000 1.0000000 +S 1 + 1 0.2067000 1.0000000 +P 3 + 1 63.4200000 0.0060440 + 2 14.6600000 0.0417990 + 3 4.4590000 0.1611430 +P 1 + 1 1.5310000 1.0000000 +P 1 + 1 0.5302000 1.0000000 +P 1 + 1 0.1750000 1.0000000 +D 1 + 1 3.7750000 1.0000000 +D 1 + 1 1.3000000 1.0000000 +D 1 + 1 0.4440000 1.0000000 +F 1 + 1 2.6660000 1.0000000 +F 1 + 1 0.8590000 1.0000000 +G 1 + 1 1.8460000 1.0000000 diff --git a/input/methods b/input/methods index 08f8a15..2aa67c1 100644 --- a/input/methods +++ b/input/methods @@ -5,7 +5,7 @@ # CCD CCSD CCSD(T) ringCCD ladderCCD F F F F F # CIS RPA RPAx ppRPA ADC - F T T F F + F F F F F # GF2 GF3 F F # G0W0 evGW qsGW diff --git a/input/molecule b/input/molecule index 779d849..4254094 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 2 7 7 0 0 # Znuc x y z - H 0. 0. 0. - H 0. 0. 1.399 + C 0. 0. 0. + O 0. 0. 2.134 diff --git a/input/molecule.xyz b/input/molecule.xyz index 072b7e1..fc30b09 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7403189714 + C 0.0000000000 0.0000000000 0.0000000000 + O 0.0000000000 0.0000000000 1.1292642494 diff --git a/input/options b/input/options index ecc7c5f..836729f 100644 --- a/input/options +++ b/input/options @@ -11,6 +11,6 @@ # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta 256 0.00001 T 5 F F T F F F F 0.000 # ACFDT: AC Kx XBS - T F F + T F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/input/weight b/input/weight index 4b80620..894420f 100644 --- a/input/weight +++ b/input/weight @@ -1,46 +1,98 @@ -1 10 -S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 +1 15 +S 9 + 1 33980.0000000 0.0000910 + 2 5089.0000000 0.0007040 + 3 1157.0000000 0.0036930 + 4 326.6000000 0.0153600 + 5 106.1000000 0.0529290 + 6 38.1100000 0.1470430 + 7 14.7500000 0.3056310 + 8 6.0350000 0.3993450 + 9 2.5300000 0.2170510 +S 9 + 1 33980.0000000 -0.0000190 + 2 5089.0000000 -0.0001510 + 3 1157.0000000 -0.0007850 + 4 326.6000000 -0.0033240 + 5 106.1000000 -0.0115120 + 6 38.1100000 -0.0341600 + 7 14.7500000 -0.0771730 + 8 6.0350000 -0.1414930 + 9 2.5300000 -0.1180190 S 1 - 1 0.7977000 1.0000000 + 1 0.7355000 1.0000000 S 1 - 1 0.2581000 1.0000000 + 1 0.2905000 1.0000000 S 1 - 1 0.0898900 1.0000000 + 1 0.1111000 1.0000000 +P 3 + 1 34.5100000 0.0053780 + 2 7.9150000 0.0361320 + 3 2.3680000 0.1424930 P 1 - 1 2.2920000 1.0000000 + 1 0.8132000 1.0000000 P 1 - 1 0.8380000 1.0000000 + 1 0.2890000 1.0000000 P 1 - 1 0.2920000 1.0000000 + 1 0.1007000 1.0000000 D 1 - 1 2.0620000 1.0000000 + 1 1.8480000 1.0000000 D 1 - 1 0.6620000 1.0000000 + 1 0.6490000 1.0000000 +D 1 + 1 0.2280000 1.0000000 F 1 - 1 1.3970000 1.0000000 -2 10 -S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 -S 1 - 1 0.7977000 1.0000000 -S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 -P 1 - 1 2.2920000 1.0000000 -P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 -D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 + 1 1.4190000 1.0000000 F 1 - 1 1.3970000 1.0000000 + 1 0.4850000 1.0000000 +G 1 + 1 1.0110000 1.0000000 +2 15 +S 9 + 1 61420.0000000 0.0000900 + 2 9199.0000000 0.0006980 + 3 2091.0000000 0.0036640 + 4 590.9000000 0.0152180 + 5 192.3000000 0.0524230 + 6 69.3200000 0.1459210 + 7 26.9700000 0.3052580 + 8 11.1000000 0.3985080 + 9 4.6820000 0.2169800 +S 9 + 1 61420.0000000 -0.0000200 + 2 9199.0000000 -0.0001590 + 3 2091.0000000 -0.0008290 + 4 590.9000000 -0.0035080 + 5 192.3000000 -0.0121560 + 6 69.3200000 -0.0362610 + 7 26.9700000 -0.0829920 + 8 11.1000000 -0.1520900 + 9 4.6820000 -0.1153310 +S 1 + 1 1.4280000 1.0000000 +S 1 + 1 0.5547000 1.0000000 +S 1 + 1 0.2067000 1.0000000 +P 3 + 1 63.4200000 0.0060440 + 2 14.6600000 0.0417990 + 3 4.4590000 0.1611430 +P 1 + 1 1.5310000 1.0000000 +P 1 + 1 0.5302000 1.0000000 +P 1 + 1 0.1750000 1.0000000 +D 1 + 1 3.7750000 1.0000000 +D 1 + 1 1.3000000 1.0000000 +D 1 + 1 0.4440000 1.0000000 +F 1 + 1 2.6660000 1.0000000 +F 1 + 1 0.8590000 1.0000000 +G 1 + 1 1.8460000 1.0000000 diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index 60a133f..8690eda 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -39,9 +39,9 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & ispin = 1 EcBSE(ispin) = 0d0 - call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) +! call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & +! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) +! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 75d35ff..58f35bb 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -45,6 +45,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:) + double precision,allocatable :: eGWlin(:) + ! Output variables double precision :: eGW(nBas) @@ -74,7 +76,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif ! Memory allocation allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & - rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin)) + rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin),eGWlin(nBas)) ! Compute linear response @@ -97,11 +99,17 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif ! Solve the quasi-particle equation - eGW(:) = eHF(:) + Z(:)*SigC(:) + eGWlin(:) = eHF(:) + Z(:)*SigC(:) ! Find all the roots of the QP equation if necessary -! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) +! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) + +! Find graphical solution of the QP equation + +! call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) + + eGW(:) = eGWlin(:) ! Dump results diff --git a/src/QuAcK/QP_graph.f90 b/src/QuAcK/QP_graph.f90 new file mode 100644 index 0000000..bc6ae48 --- /dev/null +++ b/src/QuAcK/QP_graph.f90 @@ -0,0 +1,82 @@ +subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin,eGW) + +! Compute the graphical solution of the QP equation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: eta + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + + double precision,intent(in) :: eGWlin(nBas) + +! Local variables + + integer :: p + integer :: nIt + integer,parameter :: maxIt = 64 + double precision,parameter :: thresh = 1d-6 + double precision,external :: SigmaC,dSigmaC + double precision :: sig,dsig + double precision :: f,df + double precision :: w + +! Output variables + + double precision,intent(out) :: eGW(nBas) + +! Run Newton's algorithm to find the root + + do p=nC+1,nBas-nR + + write(*,*) '-----------------' + write(*,'(A10,I3)') 'Orbital ',p + write(*,*) '-----------------' + + w = eGWlin(p) + nIt = 0 + f = 1d0 + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f + + do while (abs(f) > thresh .and. nIt < maxIt) + + nIt = nIt + 1 + + sig = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + dsig = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + f = w - eHF(p) - sig + df = 1d0 - dsig + + w = w - f/df + + write(*,'(A3,I3,A1,1X,3F15.9)') 'It.',nIt,':',w*HaToeV,f,sig + + + end do + + if(nIt == maxIt) then + + write(*,*) 'Newton root search has not converged!' + + else + + eGW(p) = w + + write(*,'(A32,F16.10)') 'Quasiparticle energy (eV) ',eGW(p)*HaToeV + write(*,*) + + end if + + end do + +end subroutine QP_graph diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index f4eaa92..0edae02 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -280,6 +280,11 @@ program QuAcK call cpu_time(start_AOtoMO) + write(*,*) + write(*,*) 'AO to MO transformation... Please be patient' + write(*,*) + + if(doSph) then ERI_MO_basis = ERI_AO_basis @@ -542,8 +547,8 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & + eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0