diff --git a/PyOptions.json b/PyOptions.json index 77b67a9..055b3aa 100644 --- a/PyOptions.json +++ b/PyOptions.json @@ -1,10 +1,10 @@ { "Scan": { - "Start":0.8, - "Stop":1.2, - "Step":0.01 + "Start":1.8, + "Stop":1.9, + "Step":0.1 }, - "Basis":"6-31G", + "Basis":"VDZ", "Outputs": { "DataOutput": { "Enabled":true, diff --git a/PyOptions.template.json b/PyOptions.template.json index 77b67a9..055b3aa 100644 --- a/PyOptions.template.json +++ b/PyOptions.template.json @@ -1,10 +1,10 @@ { "Scan": { - "Start":0.8, - "Stop":1.2, - "Step":0.01 + "Start":1.8, + "Stop":1.9, + "Step":0.1 }, - "Basis":"6-31G", + "Basis":"VDZ", "Outputs": { "DataOutput": { "Enabled":true, diff --git a/examples/molecule.N2 b/examples/molecule.N2 index ed08305..0919603 100644 --- a/examples/molecule.N2 +++ b/examples/molecule.N2 @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 7 7 0 0 + 2 7 7 4 0 # Znuc x y z N 0. 0. 0. - N 0. 0. 2.7 + N 0. 0. 2.2 diff --git a/input/basis b/input/basis index c26740e..120ad98 100644 --- a/input/basis +++ b/input/basis @@ -1,38 +1,58 @@ -1 6 +1 6 S 8 1.00 - 1469.0000000 0.0007660 - 220.5000000 0.0058920 - 50.2600000 0.0296710 - 14.2400000 0.1091800 - 4.5810000 0.2827890 - 1.5800000 0.4531230 - 0.5640000 0.2747740 - 0.0734500 0.0097510 + 9046.0000000 0.0007000 + 1357.0000000 0.0053890 + 309.3000000 0.0274060 + 87.7300000 0.1032070 + 28.5600000 0.2787230 + 10.2100000 0.4485400 + 3.8380000 0.2782380 + 0.7466000 0.0154400 S 8 1.00 - 1469.0000000 -0.0001200 - 220.5000000 -0.0009230 - 50.2600000 -0.0046890 - 14.2400000 -0.0176820 - 4.5810000 -0.0489020 - 1.5800000 -0.0960090 - 0.5640000 -0.1363800 - 0.0734500 0.5751020 + 9046.0000000 -0.0001530 + 1357.0000000 -0.0012080 + 309.3000000 -0.0059920 + 87.7300000 -0.0245440 + 28.5600000 -0.0674590 + 10.2100000 -0.1580780 + 3.8380000 -0.1218310 + 0.7466000 0.5490030 S 1 1.00 - 0.0280500 1.0000000 + 0.2248000 1.0000000 P 3 1.00 - 1.5340000 0.0227840 - 0.2749000 0.1391070 - 0.0736200 0.5003750 + 13.5500000 0.0399190 + 2.9170000 0.2171690 + 0.7973000 0.5103190 P 1 1.00 - 0.0240300 1.0000000 + 0.2185000 1.0000000 D 1 1.00 - 0.1239000 1.0000000 -2 3 -S 3 1.00 - 13.0100000 0.0196850 - 1.9620000 0.1379770 - 0.4446000 0.4781480 + 0.8170000 1.0000000 +2 6 +S 8 1.00 + 9046.0000000 0.0007000 + 1357.0000000 0.0053890 + 309.3000000 0.0274060 + 87.7300000 0.1032070 + 28.5600000 0.2787230 + 10.2100000 0.4485400 + 3.8380000 0.2782380 + 0.7466000 0.0154400 +S 8 1.00 + 9046.0000000 -0.0001530 + 1357.0000000 -0.0012080 + 309.3000000 -0.0059920 + 87.7300000 -0.0245440 + 28.5600000 -0.0674590 + 10.2100000 -0.1580780 + 3.8380000 -0.1218310 + 0.7466000 0.5490030 S 1 1.00 - 0.1220000 1.0000000 + 0.2248000 1.0000000 +P 3 1.00 + 13.5500000 0.0399190 + 2.9170000 0.2171690 + 0.7973000 0.5103190 P 1 1.00 - 0.7270000 1.0000000 + 0.2185000 1.0000000 +D 1 1.00 + 0.8170000 1.0000000 diff --git a/input/methods b/input/methods index d324776..40fd2a7 100644 --- a/input/methods +++ b/input/methods @@ -9,6 +9,6 @@ # GF2 GF3 F F # G0W0 evGW qsGW - T F T + T T F # MCMP2 F diff --git a/input/molecule b/input/molecule index f276bf5..0919603 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 2 2 0 0 + 2 7 7 4 0 # Znuc x y z - Li 0. 0. 0.000 - H 0. 0. 3.015 + N 0. 0. 0. + N 0. 0. 2.2 diff --git a/input/options b/input/options index 0373c08..dea18f7 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # RHF: maxSCF thresh DIIS n_diis guess_type ortho_type - 64 0.0000001 T 2 1 1 + 64 0.0000001 T 5 2 1 # MP: # CC: maxSCF thresh DIIS n_diis @@ -8,7 +8,7 @@ T T # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 10 3 -# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize - 128 0.00001 T 3 F F T F F F F +# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta + 256 0.00001 T 5 F F T F F F F 0.1 # 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 c26740e..120ad98 100644 --- a/input/weight +++ b/input/weight @@ -1,38 +1,58 @@ -1 6 +1 6 S 8 1.00 - 1469.0000000 0.0007660 - 220.5000000 0.0058920 - 50.2600000 0.0296710 - 14.2400000 0.1091800 - 4.5810000 0.2827890 - 1.5800000 0.4531230 - 0.5640000 0.2747740 - 0.0734500 0.0097510 + 9046.0000000 0.0007000 + 1357.0000000 0.0053890 + 309.3000000 0.0274060 + 87.7300000 0.1032070 + 28.5600000 0.2787230 + 10.2100000 0.4485400 + 3.8380000 0.2782380 + 0.7466000 0.0154400 S 8 1.00 - 1469.0000000 -0.0001200 - 220.5000000 -0.0009230 - 50.2600000 -0.0046890 - 14.2400000 -0.0176820 - 4.5810000 -0.0489020 - 1.5800000 -0.0960090 - 0.5640000 -0.1363800 - 0.0734500 0.5751020 + 9046.0000000 -0.0001530 + 1357.0000000 -0.0012080 + 309.3000000 -0.0059920 + 87.7300000 -0.0245440 + 28.5600000 -0.0674590 + 10.2100000 -0.1580780 + 3.8380000 -0.1218310 + 0.7466000 0.5490030 S 1 1.00 - 0.0280500 1.0000000 + 0.2248000 1.0000000 P 3 1.00 - 1.5340000 0.0227840 - 0.2749000 0.1391070 - 0.0736200 0.5003750 + 13.5500000 0.0399190 + 2.9170000 0.2171690 + 0.7973000 0.5103190 P 1 1.00 - 0.0240300 1.0000000 + 0.2185000 1.0000000 D 1 1.00 - 0.1239000 1.0000000 -2 3 -S 3 1.00 - 13.0100000 0.0196850 - 1.9620000 0.1379770 - 0.4446000 0.4781480 + 0.8170000 1.0000000 +2 6 +S 8 1.00 + 9046.0000000 0.0007000 + 1357.0000000 0.0053890 + 309.3000000 0.0274060 + 87.7300000 0.1032070 + 28.5600000 0.2787230 + 10.2100000 0.4485400 + 3.8380000 0.2782380 + 0.7466000 0.0154400 +S 8 1.00 + 9046.0000000 -0.0001530 + 1357.0000000 -0.0012080 + 309.3000000 -0.0059920 + 87.7300000 -0.0245440 + 28.5600000 -0.0674590 + 10.2100000 -0.1580780 + 3.8380000 -0.1218310 + 0.7466000 0.5490030 S 1 1.00 - 0.1220000 1.0000000 + 0.2248000 1.0000000 +P 3 1.00 + 13.5500000 0.0399190 + 2.9170000 0.2171690 + 0.7973000 0.5103190 P 1 1.00 - 0.7270000 1.0000000 + 0.2185000 1.0000000 +D 1 1.00 + 0.8170000 1.0000000 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 4a80e97..d1bcc76 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,4 +1,4 @@ -subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & +subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) ! Perform G0W0 calculation @@ -14,6 +14,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & logical,intent(in) :: TDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF @@ -44,6 +46,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ! Hello world + print*,eta + write(*,*) write(*,*)'************************************************' write(*,*)'| One-shot G0W0 calculation |' @@ -84,7 +88,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rhox(:,:,:,ispin)) - call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & + call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) ! COHSEX static approximation @@ -95,7 +99,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & else - call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call renormalization_factor(SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & + Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) endif diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 4d34525..8bbe609 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -65,7 +65,8 @@ program QuAcK double precision :: thresh_CC logical :: DIIS_CC - logical :: singlet_manifold,triplet_manifold + logical :: singlet_manifold + logical :: triplet_manifold integer :: maxSCF_GF,n_diis_GF,renormalization double precision :: thresh_GF @@ -74,6 +75,7 @@ program QuAcK integer :: maxSCF_GW,n_diis_GW double precision :: thresh_GW logical :: DIIS_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize + double precision :: eta integer :: nMC,nEq,nWalk,nPrint,iSeed double precision :: dt @@ -109,11 +111,11 @@ program QuAcK ! Read options for methods - call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & - maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - singlet_manifold,triplet_manifold, & - maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize, & + call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & + maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & + singlet_manifold,triplet_manifold, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -150,7 +152,7 @@ program QuAcK !------------------------------------------------------------------------ call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) - nS(:) = nO(:)*nV(:) + nS(:) = (nO(:) - nC(:))*nV(:) !------------------------------------------------------------------------ ! Read auxiliary basis set information @@ -455,7 +457,7 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & + call G0W0(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) @@ -472,7 +474,7 @@ program QuAcK if(doevGW) then call cpu_time(start_evGW) - call evGW(maxSCF_GW,thresh_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize, & + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,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_evGW) @@ -490,7 +492,7 @@ program QuAcK call cpu_time(start_qsGW) call qsGW(maxSCF_GW,thresh_GW,n_diis_GW, & - COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold, & + COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF) call cpu_time(end_qsGW) diff --git a/src/QuAcK/RHF.f90 b/src/QuAcK/RHF.f90 index 93b71e0..ef86d79 100644 --- a/src/QuAcK/RHF.f90 +++ b/src/QuAcK/RHF.f90 @@ -61,22 +61,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH ! Guess coefficients and eigenvalues - call mo_guess(nBas,Fp) - - if(guess_type == 1) then - - Fp = matmul(transpose(X),matmul(Hc,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) - c = matmul(X,cp) - - elseif(guess_type == 2) then - - call random_number(c) - - endif - - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P) ! ON(:) = 0d0 ! do i=1,nO diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index b75a598..16ca47f 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,4 +1,4 @@ -subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize, & +subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -22,6 +22,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize + double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) @@ -125,15 +126,15 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(G0W) then - call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & + call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call renormalization_factor(SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) else - call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eGW, & + call self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call renormalization_factor(SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:)) endif diff --git a/src/QuAcK/print_RHF.f90 b/src/QuAcK/print_RHF.f90 index 1cda0ec..eef7055 100644 --- a/src/QuAcK/print_RHF.f90 +++ b/src/QuAcK/print_RHF.f90 @@ -47,7 +47,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') 'MO coefficients' write(*,'(A50)') '---------------------------------------' -! call matout(nBas,nBas,cHF) + call matout(nBas,nBas,cHF) write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') 'MO energies' diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index b137d1b..18ac3e9 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,4 +1,4 @@ -subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold, & +subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF) ! Compute linear response @@ -19,6 +19,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani logical,intent(in) :: GW0 logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold + double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF @@ -142,19 +143,20 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Compute correlation part of the self-energy call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) + if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin)) if(G0W) then - call self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, & + call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call renormalization_factor(SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) else - call self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e, & + call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) + call renormalization_factor(SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) endif diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 5afc2a8..43d1455 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,8 +1,8 @@ -subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & - maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - singlet_manifold,triplet_manifold, & - maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize, & +subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & + maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & + singlet_manifold,triplet_manifold, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,renormalization, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -43,6 +43,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: G0W logical,intent(out) :: GW0 logical,intent(out) :: linearize + double precision,intent(out) :: eta integer,intent(out) :: nMC integer,intent(out) :: nEq @@ -133,10 +134,11 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t G0W = .false. GW0 = .false. linearize = .false. + eta = 0d0 read(1,*) read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, & - answer3,answer4,answer5,answer6,answer7,answer8 + answer3,answer4,answer5,answer6,answer7,answer8,eta if(answer1 == 'T') DIIS_GW = .true. if(answer2 == 'T') COHSEX = .true. diff --git a/src/QuAcK/renormalization_factor.f90 b/src/QuAcK/renormalization_factor.f90 index b4fff31..6a3e9c1 100644 --- a/src/QuAcK/renormalization_factor.f90 +++ b/src/QuAcK/renormalization_factor.f90 @@ -1,4 +1,4 @@ -subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) +subroutine renormalization_factor(SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) ! Compute renormalization factor for GW @@ -8,27 +8,25 @@ subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) ! Input variables logical,intent(in) :: SOSEX + double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: e(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS) + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables integer :: i,j,a,b,x,jb double precision :: eps - double precision,allocatable :: SigC(:),dSigC(:),d2SigC(:) - double precision,external :: Z_dcgw ! Output variables double precision,intent(out) :: Z(nBas) -! Allocate +! Initialize - allocate(SigC(nBas),dSigC(nBas),d2SigC(nBas)) - - SigC(:) = 0d0 - dSigC(:) = 0d0 - d2SigC(:) = 0d0 + Z(:) = 0d0 ! Occupied part of the correlation self-energy @@ -38,11 +36,8 @@ subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) -! Z(x) = Z(x) + 2d0*Z_dcgw(eps,rho(x,i,jb)) -! SigC(x) = SigC(x) + 2d0*rho(x,i,jb)**2/eps - dSigC(x) = dSigC(x) - 2d0*rho(x,i,jb)**2/eps**2 -! d2SigC(x) = d2SigC(x) + 4d0*rho(x,i,jb)**2/eps**3 + eps = e(x) - e(i) + Omega(jb) + Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2 enddo enddo enddo @@ -56,11 +51,8 @@ subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) -! Z(x) = Z(x) + 2d0*Z_dcgw(eps,rho(x,a,jb)) -! SigC(x) = SigC(x) + 2d0*rho(x,a,jb)**2/eps - dSigC(x) = dSigC(x) - 2d0*rho(x,a,jb)**2/eps**2 -! d2SigC(x) = d2SigC(x) + 4d0*rho(x,a,jb)**2/eps**3 + eps = e(x) - e(a) - Omega(jb) + Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2 enddo enddo enddo @@ -78,8 +70,8 @@ subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) - dSigC(x) = dSigC(x) - (rho(x,i,jb)/eps)*(rhox(x,i,jb)/eps) + eps = e(x) - e(i) + Omega(jb) + Z(x) = Z(x) - (rho(x,i,jb)/eps)*(rhox(x,i,jb)/eps) enddo enddo enddo @@ -93,8 +85,8 @@ subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) - dSigC(x) = dSigC(x) - (rho(x,a,jb)/eps)*(rhox(x,a,jb)/eps) + eps = e(x) - e(a) - Omega(jb) + Z(x) = Z(x) - (rho(x,a,jb)/eps)*(rhox(x,a,jb)/eps) enddo enddo enddo @@ -104,9 +96,6 @@ subroutine renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,Z) ! Compute renormalization factor from derivative of SigC - Z(:) = 1d0/(1d0-dSigC(:)) - -! Z(:) = 1d0 - dSigC(:) + sqrt( (1d0 - dSigC(:))**2 - 2d0*SigC(:)*d2SigC(:) ) -! Z(:) = Z(:)/(SigC(:)*d2SigC(:)) + Z(:) = 1d0/(1d0 - Z(:)) end subroutine renormalization_factor diff --git a/src/QuAcK/self_energy_correlation.f90 b/src/QuAcK/self_energy_correlation.f90 index 8d42b70..0f29b3a 100644 --- a/src/QuAcK/self_energy_correlation.f90 +++ b/src/QuAcK/self_energy_correlation.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) +subroutine self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) ! Compute correlation part of the self-energy @@ -7,14 +7,19 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, ! Input variables - logical,intent(in) :: COHSEX,SOSEX + logical,intent(in) :: COHSEX + logical,intent(in) :: SOSEX + double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: e(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS) + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables integer :: i,j,a,b,p,x,y,jb - double precision :: eps,eta + double precision :: eps ! Output variables @@ -25,11 +30,6 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, SigC = 0d0 -! Infinitesimal - - eta = 0.0d0 -! eta = 0.001d0 - ! COHSEX static approximation if(COHSEX) then @@ -118,8 +118,8 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) - SigC(x,y) = SigC(x,y) - rho(x,i,jb)*rhox(y,i,jb)/eps + eps = e(x) - e(i) + Omega(jb) + SigC(x,y) = SigC(x,y) - rho(x,i,jb)*rhox(y,i,jb)*eps/(eps**2 + eta**2) enddo enddo enddo @@ -135,8 +135,8 @@ subroutine self_energy_correlation(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho, do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) - SigC(x,y) = SigC(x,y) - rho(x,a,jb)*rhox(y,a,jb)/eps + eps = e(x) - e(a) - Omega(jb) + SigC(x,y) = SigC(x,y) - rho(x,a,jb)*rhox(y,a,jb)*eps/(eps**2 + eta**2) enddo enddo enddo diff --git a/src/QuAcK/self_energy_correlation_diag.f90 b/src/QuAcK/self_energy_correlation_diag.f90 index 0430057..ce96ebd 100644 --- a/src/QuAcK/self_energy_correlation_diag.f90 +++ b/src/QuAcK/self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) +subroutine self_energy_correlation_diag(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,rhox,EcGM,SigC) ! Compute diagonal of the correlation part of the self-energy @@ -7,14 +7,24 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega ! Input variables - logical,intent(in) :: COHSEX,SOSEX - integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: e(nBas),Omega(nS),rho(nBas,nBas,nS),rhox(nBas,nBas,nS) + logical,intent(in) :: COHSEX + logical,intent(in) :: SOSEX + double precision,intent(in) :: eta + 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) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: rhox(nBas,nBas,nS) ! Local variables integer :: i,j,a,b,p,x,jb - double precision :: eps,eta + double precision :: eps double precision,external :: SigC_dcgw ! Output variables @@ -26,11 +36,6 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega SigC = 0d0 -! Infinitesimal - - eta = 0d0 -! eta = 0.001d0 - ! COHSEX static approximation if(COHSEX) then @@ -116,7 +121,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do b=nO+1,nBas-nR jb = jb + 1 eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)/eps + EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) enddo enddo enddo @@ -133,7 +138,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do b=nO+1,nBas-nR jb = jb + 1 eps = e(x) - e(i) + Omega(jb) - SigC(x) = SigC(x) - rho(x,i,jb)*rhox(x,i,jb)/eps + SigC(x) = SigC(x) - rho(x,i,jb)*rhox(x,i,jb)*eps/(eps**2 + eta**2) enddo enddo enddo @@ -148,7 +153,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do b=nO+1,nBas-nR jb = jb + 1 eps = e(x) - e(a) - Omega(jb) - SigC(x) = SigC(x) - rho(x,a,jb)*rhox(x,a,jb)/eps + SigC(x) = SigC(x) - rho(x,a,jb)*rhox(x,a,jb)*eps/(eps**2 + eta**2) enddo enddo enddo @@ -163,7 +168,7 @@ subroutine self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,e,Omega do b=nO+1,nBas-nR jb = jb + 1 eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)/eps + EcGM = EcGM + rho(a,i,jb)*rhox(a,i,jb)*eps/(eps**2 + eta**2) enddo enddo enddo