diff --git a/GoQCaml b/GoQCaml index 5ab2190..57b0d2c 100755 --- a/GoQCaml +++ b/GoQCaml @@ -1,5 +1,5 @@ #! /bin/bash cd int -../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -###../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -m 0.5 +###../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz +../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -m 0.5 diff --git a/examples/molecule.F2 b/examples/molecule.F2 index 12bba91..9c99be9 100644 --- a/examples/molecule.F2 +++ b/examples/molecule.F2 @@ -2,4 +2,4 @@ 2 9 9 0 0 # Znuc x y z F 0. 0. 0. - F 0. 0. 3.3 + F 0. 0. 2 diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 81c624a..02328b1 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 1.4 + H 0. 0. 2.5 diff --git a/input/basis b/input/basis index fb05e68..6f3d2a9 100644 --- a/input/basis +++ b/input/basis @@ -1,18 +1,39 @@ -1 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 +1 10 +S 8 + 1 24350.0000000 0.0005020 + 2 3650.0000000 0.0038810 + 3 829.6000000 0.0199970 + 4 234.0000000 0.0784180 + 5 75.6100000 0.2296760 + 6 26.7300000 0.4327220 + 7 9.9270000 0.3506420 + 8 1.1020000 -0.0076450 +S 8 + 1 24350.0000000 -0.0001180 + 2 3650.0000000 -0.0009150 + 3 829.6000000 -0.0047370 + 4 234.0000000 -0.0192330 + 5 75.6100000 -0.0603690 + 6 26.7300000 -0.1425080 + 7 9.9270000 -0.1777100 + 8 1.1020000 0.6058360 S 1 - 1 0.1220000 1.0000000 -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 + 1 2.8360000 1.0000000 S 1 - 1 0.1220000 1.0000000 + 1 0.3782000 1.0000000 +P 3 + 1 54.7000000 0.0171510 + 2 12.4300000 0.1076560 + 3 3.6790000 0.3216810 P 1 - 1 0.7270000 1.0000000 + 1 1.1430000 1.0000000 +P 1 + 1 0.3300000 1.0000000 +D 1 + 1 4.0140000 1.0000000 +D 1 + 1 1.0960000 1.0000000 +F 1 + 1 2.5440000 1.0000000 + + diff --git a/input/methods b/input/methods index 949e2e0..bd52455 100644 --- a/input/methods +++ b/input/methods @@ -11,8 +11,8 @@ # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - F T F + T F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F diff --git a/input/molecule b/input/molecule index 81c624a..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - H 0. 0. 0. - H 0. 0. 1.4 + Ne 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..1c70680 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + Ne 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index f2ec2c1..ee86c35 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # CIS/TDHF/BSE: singlet triplet - T T + T F # GF: maxSCF thresh DIIS n_diis lin renorm 256 0.00001 T 5 T 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta diff --git a/input/weight b/input/weight index fb05e68..6f3d2a9 100644 --- a/input/weight +++ b/input/weight @@ -1,18 +1,39 @@ -1 3 -S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 +1 10 +S 8 + 1 24350.0000000 0.0005020 + 2 3650.0000000 0.0038810 + 3 829.6000000 0.0199970 + 4 234.0000000 0.0784180 + 5 75.6100000 0.2296760 + 6 26.7300000 0.4327220 + 7 9.9270000 0.3506420 + 8 1.1020000 -0.0076450 +S 8 + 1 24350.0000000 -0.0001180 + 2 3650.0000000 -0.0009150 + 3 829.6000000 -0.0047370 + 4 234.0000000 -0.0192330 + 5 75.6100000 -0.0603690 + 6 26.7300000 -0.1425080 + 7 9.9270000 -0.1777100 + 8 1.1020000 0.6058360 S 1 - 1 0.1220000 1.0000000 -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 + 1 2.8360000 1.0000000 S 1 - 1 0.1220000 1.0000000 + 1 0.3782000 1.0000000 +P 3 + 1 54.7000000 0.0171510 + 2 12.4300000 0.1076560 + 3 3.6790000 0.3216810 P 1 - 1 0.7270000 1.0000000 + 1 1.1430000 1.0000000 +P 1 + 1 0.3300000 1.0000000 +D 1 + 1 4.0140000 1.0000000 +D 1 + 1 1.0960000 1.0000000 +F 1 + 1 2.5440000 1.0000000 + + diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 526959d..2e19b5c 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -1,5 +1,5 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) + linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -52,8 +52,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m double precision,allocatable :: SigT(:) double precision,allocatable :: Z(:) - double precision,allocatable :: eG0T0(:) - double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) @@ -61,6 +59,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Output variables + double precision,intent(out) :: eG0T0(nBas) + ! Hello world write(*,*) @@ -88,7 +88,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), & - SigT(nBas),Z(nBas),eG0T0(nBas)) + SigT(nBas),Z(nBas)) !---------------------------------------------- ! alpha-beta block diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 0763450..8bf40b4 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -34,6 +34,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & ! Local variables + logical :: print_W = .false. integer :: ispin double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) @@ -104,7 +105,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & eGWlin(:) = eHF(:) + Z(:)*SigC(:) if(linearize) then - + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + eGW(:) = eGWlin(:) ! Find all the roots of the QP equation if necessary @@ -121,7 +125,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & ! Dump results -! call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + if(print_W) call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) ! Compute the RPA correlation energy diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index da24f7a..7c8c831 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -14,6 +14,7 @@ program QuAcK logical :: doG0W0,doevGW,doqsGW logical :: doG0T0,doevGT,doqsGT logical :: doMCMP2,doMinMCMP2 + logical :: doGTGW = .false. logical :: doBas integer :: nNuc,nBas,nBasCABS @@ -46,8 +47,10 @@ program QuAcK double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:) double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:) - double precision,allocatable :: ERI_AO_basis(:,:,:,:) - double precision,allocatable :: ERI_MO_basis(:,:,:,:) + double precision,allocatable :: ERI_AO(:,:,:,:) + double precision,allocatable :: ERI_MO(:,:,:,:) + double precision,allocatable :: ERI_ERF_AO(:,:,:,:) + double precision,allocatable :: ERI_ERF_MO(:,:,:,:) double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:) double precision :: start_QuAcK ,end_QuAcK ,t_QuAcK @@ -202,7 +205,7 @@ program QuAcK allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),eG0T0(nBas),PHF(nBas,nBas,nspin), & S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), & - ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas)) + ERI_AO(nBas,nBas,nBas,nBas),ERI_MO(nBas,nBas,nBas,nBas)) ! Read integrals @@ -210,12 +213,12 @@ program QuAcK if(doSph) then - call read_integrals_sph(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis) + call read_integrals_sph(nEl(:),nBas,S,T,V,Hc,ERI_AO) else call system('./GoQCaml') - call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis) + call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO) end if @@ -237,7 +240,7 @@ program QuAcK if(doRHF) then call cpu_time(start_HF) - call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO_basis,X,ENuc,ERHF,eHF,cHF,PHF) + call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,ERHF,eHF,cHF,PHF) call cpu_time(end_HF) t_HF = end_HF - start_HF @@ -253,7 +256,7 @@ program QuAcK if(doUHF) then call cpu_time(start_HF) - call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO_basis,X,ENuc,EUHF,eHF,cHF,PHF) + call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,EUHF,eHF,cHF,PHF) call cpu_time(end_HF) t_HF = end_HF - start_HF @@ -270,7 +273,7 @@ program QuAcK call cpu_time(start_MOM) call MOM(maxSCF_HF,thresh_HF,n_diis_HF, & - nBas,nO,S,T,V,Hc,ERI_AO_basis,X,ENuc,ERHF,cHF,eHF,PHF) + nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,ERHF,cHF,eHF,PHF) call cpu_time(end_MOM) t_MOM = end_MOM - start_MOM @@ -285,7 +288,7 @@ program QuAcK ! Compute Hartree Hamiltonian in the MO basis - call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H) + call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO,H) call cpu_time(start_AOtoMO) @@ -296,13 +299,13 @@ program QuAcK if(doSph) then - ERI_MO_basis = ERI_AO_basis + ERI_MO(:,:,:,:) = ERI_AO(:,:,:,:) print*,'!!! MO = AO !!!' - deallocate(ERI_AO_basis) + deallocate(ERI_AO) else - call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis) + call AOtoMO_integral_transform(nBas,cHF,ERI_AO,ERI_MO) end if @@ -319,7 +322,7 @@ program QuAcK if(doMP2) then call cpu_time(start_MP2) - call MP2(nBas,nC,nO,nV,nR,ERI_MO_basis,ENuc,ERHF,eHF,EcMP2) + call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2) call cpu_time(end_MP2) t_MP2 = end_MP2 - start_MP2 @@ -335,7 +338,7 @@ program QuAcK if(doMP3) then call cpu_time(start_MP3) - call MP3(nBas,nEl,ERI_MO_basis,eHF,ENuc,ERHF) + call MP3(nBas,nEl,ERI_MO,eHF,ENuc,ERHF) call cpu_time(end_MP3) t_MP3 = end_MP3 - start_MP3 @@ -358,8 +361,8 @@ program QuAcK ! Read integrals - call read_F12_integrals(nBas,S,ERI_AO_basis,F12,Yuk,FC) - call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,eHF,cHF) + call read_F12_integrals(nBas,S,ERI_AO,F12,Yuk,FC) + call MP2F12(nBas,nC,nO,nV,ERI_AO,F12,Yuk,FC,ERHF,eHF,cHF) call cpu_time(end_MP2F12) t_MP2F12 = end_MP2F12 - start_MP2F12 @@ -376,7 +379,7 @@ program QuAcK call cpu_time(start_CCD) call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & - ERI_MO_basis,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -395,7 +398,7 @@ program QuAcK call cpu_time(start_CCSD) call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC(1),nO(1),nV(1),nR(1), & - ERI_MO_basis,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCSD) t_CCSD = end_CCSD - start_CCSD @@ -412,7 +415,7 @@ program QuAcK call cpu_time(start_CCD) call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & - ERI_MO_basis,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -429,7 +432,7 @@ program QuAcK call cpu_time(start_CCD) call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & - ERI_MO_basis,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -446,7 +449,7 @@ program QuAcK call cpu_time(start_CCD) call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & - ERI_MO_basis,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -463,7 +466,7 @@ program QuAcK call cpu_time(start_CCD) call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), & - ERI_MO_basis,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD @@ -479,7 +482,7 @@ program QuAcK if(doCIS) then call cpu_time(start_CIS) - call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eHF) + call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF) call cpu_time(end_CIS) t_CIS = end_CIS - start_CIS @@ -496,7 +499,7 @@ program QuAcK call cpu_time(start_RPA) call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -513,7 +516,7 @@ program QuAcK call cpu_time(start_RPAx) call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPAx) t_RPAx = end_RPAx - start_RPAx @@ -530,7 +533,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_basis,eHF) + nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_ppRPA) t_ppRPA = end_ppRPA - start_ppRPA @@ -547,7 +550,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_basis) + nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO) call cpu_time(end_ADC) t_ADC = end_ADC - start_ADC @@ -563,7 +566,7 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF) + call G0F2(linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -579,7 +582,7 @@ program QuAcK if(doevGF2) then call cpu_time(start_GF2) - call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF) + call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -595,7 +598,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_basis,eHF) + call G0F3(renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) call cpu_time(end_GF3) t_GF3 = end_GF3 - start_GF3 @@ -611,7 +614,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_basis,eHF) + call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF) call cpu_time(end_GF3) t_GF3 = end_GF3 - start_GF3 @@ -631,7 +634,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & singlet_manifold,triplet_manifold,linGW,eta, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) t_G0W0 = end_G0W0 - start_G0W0 @@ -649,7 +652,7 @@ program QuAcK call cpu_time(start_evGW) call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & 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) + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) t_evGW = end_evGW - start_evGW @@ -667,7 +670,7 @@ program QuAcK call cpu_time(start_qsGW) call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,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,ERI_MO_basis,PHF,cHF,eHF) + 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) call cpu_time(end_qsGW) t_qsGW = end_qsGW - start_qsGW @@ -687,8 +690,7 @@ program QuAcK call cpu_time(start_G0T0) call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & singlet_manifold,triplet_manifold,linGW,eta, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF) - + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_G0T0) t_G0T0 = end_G0T0 - start_G0T0 @@ -705,7 +707,7 @@ program QuAcK call cpu_time(start_evGT) call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, & - eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0) + eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_evGT) t_evGT = end_evGT - start_evGT @@ -789,6 +791,64 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Range-separeted GT/GW +!------------------------------------------------------------------------ + + if(doGTGW) then + + ! Read and transform long-range two-electron integrals + + allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas)) + call read_LR(nBas,ERI_ERF_AO) + + call cpu_time(start_AOtoMO) + + write(*,*) + write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient' + write(*,*) + + call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO) + + call cpu_time(end_AOtoMO) + + deallocate(ERI_ERF_AO) + + t_AOtoMO = end_AOtoMO - start_AOtoMO + + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,*) + + ! Long-range G0W0 calculation + + call cpu_time(start_G0W0) + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & + singlet_manifold,triplet_manifold,linGW,eta, & + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) + call cpu_time(end_G0W0) + + t_G0W0 = end_G0W0 - start_G0W0 + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds' + write(*,*) + + ! Short-range G0T0 calculation + + ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:) + + call cpu_time(start_G0T0) + call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & + singlet_manifold,triplet_manifold,linGW,eta, & + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) + call cpu_time(end_G0T0) + + t_G0T0 = end_G0T0 - start_G0T0 + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' + write(*,*) + + call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV) + + end if + !------------------------------------------------------------------------ ! Basis set correction !------------------------------------------------------------------------ @@ -799,7 +859,7 @@ program QuAcK call cpu_time(start_Bas) call basis_correction(nBas,nO,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, & - ERI_MO_basis,eHF,cHF,PHF,eG0W0) + ERI_MO,eHF,cHF,PHF,eG0W0) call cpu_time(end_Bas) t_Bas = end_Bas - start_Bas diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 index b1ba603..82f00f7 100644 --- a/src/QuAcK/linear_response_C_pp.f90 +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -25,7 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) ! Define the chemical potential eF = e(nO) + e(nO+1) - eF = 0d0 +! eF = 0d0 ! Build C matrix for the singlet manifold diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 index 0a85d64..4257795 100644 --- a/src/QuAcK/linear_response_D_pp.f90 +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -25,7 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) ! Define the chemical potential eF = e(nO) + e(nO+1) - eF = 0d0 +! eF = 0d0 ! Build the D matrix for the singlet manifold diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index a709673..640d95e 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -121,10 +121,16 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) +! allocate(order(nOO+nVV)) +! call quick_sort(Omega(:),order(:),nOO+nVV) +! call matout(nOO+nVV,1,Omega(:)*HaToeV) + ! Split the various quantities in p-p and h-h parts call sort_ppRPA(ortho_eigvec,nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:)) +! call matout(32,1,(Omega1(:) - Omega1(1))*HaToeV) + ! Compute the RPA correlation energy EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index 3c45cc5..7c10627 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -56,6 +56,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER nOO = nO*(nO+1)/2 nVV = nV*(nV+1)/2 +! print*,'nOO = ',nOO +! print*,'nVV = ',nVV + ! Memory allocation allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & diff --git a/src/utils/read_LR.f90 b/src/utils/read_LR.f90 new file mode 100644 index 0000000..554b089 --- /dev/null +++ b/src/utils/read_LR.f90 @@ -0,0 +1,56 @@ +subroutine read_LR(nBas,G) + +! Read the long-range two-electron integrals from files + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + +! Local variables + + integer :: mu,nu,la,si + double precision :: ERI + double precision :: lambda + +! Output variables + + double precision,intent(out) :: G(nBas,nBas,nBas,nBas) + +! Open file with integrals + + lambda = 1d0 + + print*, 'Scaling integrals by ',lambda + + open(unit=11,file='int/ERI_lr.dat') + +! Read two-electron integrals + + G(:,:,:,:) = 0d0 + do + read(11,*,end=11) mu,nu,la,si,ERI + + ERI = lambda*ERI +! <12|34> + G(mu,nu,la,si) = ERI +! <32|14> + G(la,nu,mu,si) = ERI +! <14|32> + G(mu,si,la,nu) = ERI +! <34|12> + G(la,si,mu,nu) = ERI +! <41|23> + G(si,mu,nu,la) = ERI +! <23|41> + G(nu,la,si,mu) = ERI +! <21|43> + G(nu,mu,si,la) = ERI +! <43|21> + G(si,la,nu,mu) = ERI + enddo + 11 close(unit=11) + +end subroutine read_LR diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index 4004d1f..9ae5458 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -19,7 +19,11 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) ! Output variables - double precision,intent(out) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas) + double precision,intent(out) :: S(nBas,nBas) + double precision,intent(out) :: T(nBas,nBas) + double precision,intent(out) :: V(nBas,nBas) + double precision,intent(out) :: Hc(nBas,nBas) + double precision,intent(out) :: G(nBas,nBas,nBas,nBas) ! Open file with integrals diff --git a/src/utils/sort.f90 b/src/utils/sort.f90 new file mode 100644 index 0000000..91d4902 --- /dev/null +++ b/src/utils/sort.f90 @@ -0,0 +1,89 @@ + subroutine quick_sort(x,iorder,isize) + + implicit none + + integer,intent(in) :: isize + double precision,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + + call rec_quicksort(x,iorder,isize,1,isize,1) + +end subroutine quick_sort + +recursive subroutine rec_quicksort(x,iorder,isize,first,last,level) + + implicit none + + integer, intent(in) :: isize, first, last, level + integer,intent(inout) :: iorder(isize) + double precision, intent(inout):: x(isize) + double precision :: c, tmp + integer :: itmp + integer :: i, j + + c = x( shiftr(first+last,1) ) + i = first + j = last + do + do while (x(i) < c) + i=i+1 + end do + do while (c < x(j)) + j=j-1 + end do + if (i >= j) then + exit + endif + tmp = x(i) + x(i) = x(j) + x(j) = tmp + itmp = iorder(i) + iorder(i) = iorder(j) + iorder(j) = itmp + i=i+1 + j=j-1 + enddo + if ( ((i-first <= 10000).and.(last-j <= 10000)).or.(level<=0) ) then + if (first < i-1) then + call rec_quicksort(x, iorder, isize, first, i-1,level/2) + endif + if (j+1 < last) then + call rec_quicksort(x, iorder, isize, j+1, last,level/2) + endif + else + if (first < i-1) then + call rec_quicksort(x, iorder, isize, first, i-1,level/2) + endif + if (j+1 < last) then + call rec_quicksort(x, iorder, isize, j+1, last,level/2) + endif + endif +end subroutine rec_quicksort + +subroutine set_order(x,iorder,isize,jsize) + + implicit none + + integer :: isize,jsize + double precision :: x(isize,jsize) + double precision,allocatable :: xtmp(:,:) + integer :: iorder(*) + integer :: i,j + + allocate(xtmp(isize,jsize)) + + do i=1,isize + do j=1,jsize + xtmp(i,j) = x(i,iorder(j)) + end do + end do + + do i=1,isize + do j=1,jsize + x(i,j) = xtmp(i,j) + end do + end do + + deallocate(xtmp) + +end subroutine set_order