From ce7d0118ab6c4bafa510afdefc70053e1afb16c6 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 7 May 2019 22:55:36 +0200 Subject: [PATCH] sph ready --- GoSph | 11 ++-- extract_sph.sh | 106 +++++++++++++-------------------- input/basis | 40 +------------ input/methods | 8 +-- input/molecule | 2 +- input/options | 4 +- run_sph.sh | 10 ++-- src/QuAcK/G0W0.f90 | 27 ++++----- src/QuAcK/Makefile | 13 +++- src/QuAcK/QuAcK.f90 | 69 +++++++++++++++++---- src/QuAcK/evGW.f90 | 13 ++-- src/QuAcK/print_G0W0.f90 | 14 ++--- src/QuAcK/print_evGW.f90 | 14 ++--- src/QuAcK/print_excitation.f90 | 7 ++- src/QuAcK/print_qsGW.f90 | 4 +- src/QuAcK/qsGW.f90 | 29 ++++++--- src/utils/Makefile | 2 +- src/utils/read_integrals.f90 | 14 ++--- 18 files changed, 191 insertions(+), 196 deletions(-) diff --git a/GoSph b/GoSph index 282b774..db90c20 100755 --- a/GoSph +++ b/GoSph @@ -9,9 +9,10 @@ if [ $# = 2 ] then cp examples/molecule.Sph_"$1" input/molecule cp examples/basis.Sph.Ylm"$2" input/basis -cp ~/Integrals/QuAcK_Sph/Sph_ERI_"$2".dat int/ERI.dat -cp ~/Integrals/QuAcK_Sph/Sph_Kin_"$2".dat int/Kin.dat -cp ~/Integrals/QuAcK_Sph/Sph_Nuc_"$2".dat int/Nuc.dat -cp ~/Integrals/QuAcK_Sph/Sph_Ov_"$2".dat int/Ov.dat -./bin/QuAcK +cp ~/Integrals/QuAcK_Sph/Sph_ERI_"$2".dat ~/Integrals/QuAcK_Sph/ERI.dat +cp ~/Integrals/QuAcK_Sph/Sph_Kin_"$2".dat ~/Integrals/QuAcK_Sph/Kin.dat +cp ~/Integrals/QuAcK_Sph/Sph_Nuc_"$2".dat ~/Integrals/QuAcK_Sph/Nuc.dat +cp ~/Integrals/QuAcK_Sph/Sph_Ov_"$2".dat ~/Integrals/QuAcK_Sph/Ov.dat +./bin/QuAcK | tee sph.out + ./extract_sph.sh sph.out fi diff --git a/extract_sph.sh b/extract_sph.sh index a833748..1714728 100755 --- a/extract_sph.sh +++ b/extract_sph.sh @@ -1,86 +1,64 @@ #! /bin/bash -for i in *.out -do +INPUT=$1 echo - echo '***********************' - echo '*** ' $i ' ***' - echo '***********************' + echo '******************************************' + echo '*** Extracting information of' $INPUT ' ***' + echo '******************************************' echo echo - echo '--- HF information ---' - echo - - grep "Hartree-Fock energy" $i - grep "HF HOMO energy (eV):" $i - grep "HF LUMO energy (eV):" $i - grep "HF HOMO-LUMO gap (eV):" $i + echo '*** WFT information ***' + grep "MP2 correlation energy" $INPUT + grep "Ec(MP2) =" $INPUT + grep "Ec(CCSD) =" $INPUT + grep "Ec(CCSD(T)) =" $INPUT echo - echo '--- WFT information ---' - echo + echo '*** Gap information: HF, G0F2, GF2, G0W0 & evGW ***' + HF=`grep "HF HOMO-LUMO gap (eV):" $INPUT | cut -f2 -d":"` + G0F2=`grep "GF2 HOMO-LUMO gap (eV):" $INPUT | head -1 | cut -f2 -d":"` + GF2=`grep "GF2 HOMO-LUMO gap (eV):" $INPUT | tail -1 | cut -f2 -d":"` + G0W0=`grep "G0W0 HOMO-LUMO gap (eV):" $INPUT | cut -f2 -d":"` + evGW=`grep "evGW HOMO-LUMO gap (eV):" $INPUT | tail -1 | cut -f2 -d":"` - grep "Ec(MP2) =" $i - grep "Ec(CCSD) =" $i - grep "Ec(CCSD(T)) =" $i + echo -e "\t" $HF "\t" $G0F2 "\t" $GF2 "\t" $G0W0 "\t" $evGW echo - echo '--- CIS excitation energy (singlet & triplet) ---' - echo - grep "| 1 |" $i | head -1 | cut -f4 -d"|" - grep "| 1 |" $i | head -2 | cut -f4 -d"|" | tail -1 + echo '*** Ec@G0W0 information: RPA, GM, BSE1 & BSE3 ***' + RPA_G0W0=`grep "RPA@G0W0 correlation energy =" $INPUT| cut -f2 -d"="` + GM_G0W0=`grep "GM@G0W0 correlation energy =" $INPUT| cut -f2 -d"="` + BSE1_G0W0=`grep "BSE@G0W0 correlation energy (singlet)" $INPUT| cut -f2 -d"="` + BSE3_G0W0=`grep "BSE@G0W0 correlation energy (triplet)" $INPUT| cut -f2 -d"="` + + echo -e "\t" $RPA_G0W0 "\t" $GM_G0W0 "\t" $BSE1_G0W0 "\t" $BSE3_G0W0 echo - echo '--- TDHF excitation energy (singlet & triplet) ---' - echo - grep "| 1 |" $i | head -3 | cut -f4 -d"|" | tail -1 - grep "| 1 |" $i | head -4 | cut -f4 -d"|" | tail -1 + echo '*** Ec@evGW information: RPA, GM, BSE1 & BSE3 ***' + RPA_evGW=`grep "RPA@evGW correlation energy =" $INPUT | tail -1| cut -f2 -d"="` + GM_evGW=`grep "GM@evGW correlation energy =" $INPUT | tail -1 | cut -f2 -d"="` + BSE1_evGW=`grep "BSE@evGW correlation energy (singlet)" $INPUT | cut -f2 -d"="` + BSE3_evGW=`grep "BSE@evGW correlation energy (triplet)" $INPUT | cut -f2 -d"="` + + echo -e "\t" $RPA_evGW "\t" $GM_evGW "\t" $BSE1_evGW "\t" $BSE3_evGW echo - echo '--- GF2 information ---' - echo - - grep "GF2 HOMO energy (eV):" $i | tail -1 - grep "GF2 LUMO energy (eV):" $i | tail -1 - grep "GF2 HOMO-LUMO gap (eV):" $i | tail -1 + echo '*** CIS and TDHF excitation energy (singlet & triplet) ***' + CIS1=`grep "| 1 |" $INPUT | head -1 | cut -f4 -d"|"` + CIS3=`grep "| 1 |" $INPUT | head -2 | cut -f4 -d"|" | tail -1` + TDHF1=`grep "| 1 |" $INPUT | head -3 | cut -f4 -d"|" | tail -1` + TDHF3=`grep "| 1 |" $INPUT | head -4 | cut -f4 -d"|" | tail -1` + echo -e "\t" $CIS1 "\t" $CIS3 "\t" $TDHF1 "\t" $TDHF3 echo - echo '--- G0W0 information ---' - echo - - grep "G0W0 HOMO energy (eV):" $i - grep "G0W0 LUMO energy (eV):" $i - grep "G0W0 HOMO-LUMO gap (eV):" $i - grep "G0W0 RPA total energy =" $i - grep "G0W0 GM total energy =" $i - - echo - echo '--- BSE@G0W0 excitation energy (singlet & triplet) ---' - echo - grep "| 1 |" $i | head -6 | cut -f4 -d"|" | tail -1 - grep "| 1 |" $i | head -7 | cut -f4 -d"|" | tail -1 - - echo - echo '--- evGW information ---' - echo - - grep "evGW HOMO energy (eV):" $i | tail -1 - grep "evGW LUMO energy (eV):" $i | tail -1 - grep "evGW HOMO-LUMO gap (eV):" $i | tail -1 - grep "evGW RPA total energy =" $i | tail -1 - grep "evGW GM total energy =" $i | tail -1 - - echo - echo '--- BSE@evGW excitation energy (singlet & triplet) ---' - echo - grep "| 1 |" $i | tail -2 | head -1 | cut -f4 -d"|" - grep "| 1 |" $i | tail -1 | cut -f4 -d"|" + echo '*** BSE@G0W0 and BSE@evGW excitation energy (singlet & triplet) ***' + G0W01=`grep "| 1 |" $INPUT | head -6 | cut -f4 -d"|" | tail -1` + G0W03=`grep "| 1 |" $INPUT | head -7 | cut -f4 -d"|" | tail -1` + evGW1=`grep "| 1 |" $INPUT | tail -2 | head -1 | cut -f4 -d"|"` + evGW3=`grep "| 1 |" $INPUT | tail -1 | cut -f4 -d"|"` + echo -e "\t" $G0W01 "\t" $G0W03 "\t" $evGW1 "\t" $evGW3 echo echo '*** DONE ***' - echo - -done diff --git a/input/basis b/input/basis index 18db43a..8a72e6e 100644 --- a/input/basis +++ b/input/basis @@ -1,42 +1,4 @@ -1 100 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 -S 1 1.00 - 1.0000000 1.0000000 +1 81 S 1 1.00 1.0000000 1.0000000 S 1 1.00 diff --git a/input/methods b/input/methods index 341e057..ca18179 100644 --- a/input/methods +++ b/input/methods @@ -1,14 +1,14 @@ # RHF UHF MOM T F F # MP2 MP3 MP2-F12 - F F F + T F F # CCD CCSD CCSD(T) F F F # CIS TDHF ADC - F F F + T T F # GF2 GF3 - F F + T F # G0W0 evGW qsGW - T F F + T T F # MCMP2 F diff --git a/input/molecule b/input/molecule index 321523e..50795c8 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 9 9 0 0 + 1 16 16 0 0 # Znuc x y z X 0.0 0.0 0.0 diff --git a/input/options b/input/options index 5095ea0..0113ebb 100644 --- a/input/options +++ b/input/options @@ -7,8 +7,8 @@ # CIS/TDHF: singlet triplet T T # GF: maxSCF thresh DIIS n_diis renormalization - 64 0.00001 T 5 3 + 64 0.00001 T 10 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize - 64 0.00001 T 3 F F T F F F F + 64 0.00001 T 10 F F T F F F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/run_sph.sh b/run_sph.sh index bfe1c48..fd1cc87 100755 --- a/run_sph.sh +++ b/run_sph.sh @@ -1,8 +1,8 @@ #! /bin/bash -Lmin=0 -Lmax=0 -Mmax=9 +Lmin=1 +Lmax=1 +Mmax=10 rs=$1 if [ $# != 1 ] @@ -30,8 +30,8 @@ else nb=$(bc -l <<< "(($M+1)*($M+1))") echo "Number of basis functions = " $nb echo -e "# rs \n" $rs > input/sph - ./GoSph $ne $M > Sph_${ne}_${nb}.out - grep "Total CPU time for QuAcK =" Sph_${ne}_${nb}.out + ./GoSph $ne $M > out/Sph_${ne}_${nb}.out + grep "Total CPU time for QuAcK =" out/Sph_${ne}_${nb}.out done diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 25770ed..f7b4394 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,5 +1,5 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) ! Perform G0W0 calculation @@ -20,8 +20,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO_basis(nBas,nBas,nBas,nBas) + double precision,intent(in) :: H(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables @@ -30,7 +30,6 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcGM - double precision,allocatable :: H(:,:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) double precision,allocatable :: Omega(:,:) @@ -65,23 +64,19 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ! Memory allocation - allocate(H(nBas,nBas),SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin), & + allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin), & rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin)) -! Compute Hartree Hamiltonian in the MO basis - - call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H) - ! Compute linear response - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, & + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) ! Compute correlation part of the self-energy - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin)) + 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, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) @@ -122,7 +117,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ispin = 1 EcBSE(ispin) = 0d0 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, & + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI, & rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) @@ -135,11 +130,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & ispin = 2 EcBSE(ispin) = 0d0 - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, & + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, & + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI, & rho(:,:,:,1),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) diff --git a/src/QuAcK/Makefile b/src/QuAcK/Makefile index 2c54de3..5183c73 100644 --- a/src/QuAcK/Makefile +++ b/src/QuAcK/Makefile @@ -3,11 +3,15 @@ BDIR =../../bin ODIR = obj OODIR = ../IntPak/obj SDIR =. -FC = gfortran -I$(IDIR) +FC = gfortran -I$(IDIR) -Wall -g -Wno-unused -Wno-unused-dummy-argument ifeq ($(DEBUG),1) -FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant + FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant else -FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2 + FFLAGS = -O3 +endif + +ifeq ($(PROFIL),1) + FC += -p -fno-inline endif LIBS = ~/Dropbox/quack/lib/*.a @@ -32,5 +36,8 @@ debug: DEBUG=1 make $(BDIR)/QuAcK #DEBUG=1 make clean $(BDIR)/QuAcK +profil: + PROFIL=1 make $(BDIR)/QuAcK + clean: rm -f $(ODIR)/*.o $(BDIR)/QuAcK $(BDIR)/debug diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 18e82ad..310f82c 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -3,6 +3,7 @@ program QuAcK implicit none include 'parameters.h' + logical :: doSph logical :: doRHF,doUHF,doMOM logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT @@ -27,13 +28,15 @@ program QuAcK integer :: TrialType double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:) - double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:) + double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:) double precision,allocatable :: ERI_AO_basis(:,:,:,:),ERI_MO_basis(:,:,:,:) double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:) double precision :: start_QuAcK ,end_QuAcK ,t_QuAcK + double precision :: start_int ,end_int ,t_int double precision :: start_HF ,end_HF ,t_HF double precision :: start_MOM ,end_MOM ,t_MOM + double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO double precision :: start_CCD ,end_CCD ,t_CCD double precision :: start_CCSD ,end_CCSD ,t_CCSD double precision :: start_CIS ,end_CIS ,t_CIS @@ -85,6 +88,10 @@ program QuAcK write(*,*) '******************************************************************************************' write(*,*) +! Spherium calculation? + + doSph = .true. + call cpu_time(start_QuAcK) ! Which calculations do you want to do? @@ -162,13 +169,30 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),PHF(nBas,nBas,nspin), & - S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), & + allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(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)) ! Read integrals - call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis) + call cpu_time(start_int) + + if(doSph) then + + call read_integrals_sph(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis) + + else + + call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis) + + end if + + call cpu_time(end_int) + + t_int = end_int - start_int + write(*,*) + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for reading integrals = ',t_int,' seconds' + write(*,*) ! Compute orthogonalization matrix @@ -185,7 +209,7 @@ program QuAcK call cpu_time(end_HF) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for UHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds' write(*,*) endif @@ -227,9 +251,29 @@ program QuAcK ! AO to MO integral transform for post-HF methods !------------------------------------------------------------------------ -! call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis) - ERI_MO_basis = ERI_AO_basis - print*,'!!! MO = AO !!!' +! Compute Hartree Hamiltonian in the MO basis + + call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H) + + call cpu_time(start_AOtoMO) + + if(doSph) then + + ERI_MO_basis = ERI_AO_basis + print*,'!!! MO = AO !!!' + deallocate(ERI_AO_basis) + + else + + call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis) + + end if + + call cpu_time(end_AOtoMO) + + t_AOtoMO = end_AOtoMO - start_AOtoMO + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,*) !------------------------------------------------------------------------ ! Compute MP2 energy @@ -326,8 +370,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_basis,eHF) call cpu_time(end_CIS) t_CIS = end_CIS - start_CIS @@ -410,7 +453,7 @@ program QuAcK call cpu_time(start_G0W0) call G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + 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 @@ -427,7 +470,7 @@ program QuAcK 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, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + 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) t_evGW = end_evGW - start_evGW @@ -445,7 +488,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, & - nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF) + 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) t_qsGW = end_qsGW - start_qsGW diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index e7443c0..7386484 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,5 +1,5 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -28,7 +28,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,intent(in) :: PHF(nBas,nBas) double precision,intent(in) :: eG0W0(nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas) + double precision,intent(in) :: H(nBas,nBas) double precision,intent(in) :: ERI_MO_basis(nBas,nBas,nBas,nBas) ! Local variables @@ -49,7 +49,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,allocatable :: eGW(:) double precision,allocatable :: eOld(:) double precision,allocatable :: Z(:) - double precision,allocatable :: H(:,:) double precision,allocatable :: SigC(:) double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) @@ -80,8 +79,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Memory allocation - allocate(eGW(nBas),eOld(nBas),Z(nBas),H(nBas,nBas),SigC(nBas),Omega(nS,nspin), & - XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & + allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), & + XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Initialization @@ -96,10 +95,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani eOld(:) = eGW(:) Z(:) = 1d0 -! Compute Hartree Hamiltonian in the MO basis - - call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H) - !------------------------------------------------------------------------ ! Main loop !------------------------------------------------------------------------ diff --git a/src/QuAcK/print_G0W0.f90 b/src/QuAcK/print_G0W0.f90 index d16769a..782643b 100644 --- a/src/QuAcK/print_G0W0.f90 +++ b/src/QuAcK/print_G0W0.f90 @@ -36,14 +36,14 @@ subroutine print_G0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) enddo write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A27,F15.6)') 'G0W0 HOMO energy (eV):',eGW(HOMO)*HaToeV - write(*,'(2X,A27,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV - write(*,'(2X,A27,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO energy (eV):',eGW(HOMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 LUMO energy (eV):',eGW(LUMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A27,F15.6)') 'G0W0 RPA total energy =',ENuc + EHF + EcRPA - write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA - write(*,'(2X,A27,F15.6)') 'G0W0 GM total energy =',ENuc + EHF + EcGM - write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@G0W0 correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@G0W0 total energy =',ENuc + EHF + EcGM + write(*,'(2X,A30,F15.6)') 'GM@G0W0 correlation energy =',EcGM write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/print_evGW.f90 b/src/QuAcK/print_evGW.f90 index eaed01f..fe61210 100644 --- a/src/QuAcK/print_evGW.f90 +++ b/src/QuAcK/print_evGW.f90 @@ -43,14 +43,14 @@ subroutine print_evGW(nBas,nO,nSCF,Conv,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA,EcGM) write(*,'(2X,A10,I3)') 'Iteration ',nSCF write(*,'(2X,A14,F15.5)')'Convergence = ',Conv write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A27,F15.6)') 'evGW HOMO energy (eV):',eGW(HOMO)*HaToeV - write(*,'(2X,A27,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV - write(*,'(2X,A27,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW HOMO energy (eV):',eGW(HOMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW LUMO energy (eV):',eGW(LUMO)*HaToeV + write(*,'(2X,A30,F15.6)') 'evGW HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A27,F15.6)') 'evGW RPA total energy =',ENuc + EHF + EcRPA - write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA - write(*,'(2X,A27,F15.6)') 'evGW GM total energy =',ENuc + EHF + EcGM - write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM + write(*,'(2X,A30,F15.6)') 'RPA@evGW total energy =',ENuc + EHF + EcRPA + write(*,'(2X,A30,F15.6)') 'RPA@evGW correlation energy =',EcRPA + write(*,'(2X,A30,F15.6)') 'GM@evGW total energy =',ENuc + EHF + EcGM + write(*,'(2X,A30,F15.6)') 'GM@evGW correlation energy =',EcGM write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index dbf90ef..889ed6e 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -5,11 +5,16 @@ subroutine print_excitation(method,ispin,nS,Omega) implicit none include 'parameters.h' +! Input variables + character*5,intent(in) :: method integer,intent(in) :: ispin,nS double precision,intent(in) :: Omega(nS) +! Local variables + character*7 :: spin_manifold + integer,parameter :: maxS = 32 integer :: ia if(ispin == 1) spin_manifold = 'singlet' @@ -23,7 +28,7 @@ subroutine print_excitation(method,ispin,nS,Omega) '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|' write(*,*)'-------------------------------------------------------------' - do ia=1,nS + do ia=1,min(nS,maxS) write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & '|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|' enddo diff --git a/src/QuAcK/print_qsGW.f90 b/src/QuAcK/print_qsGW.f90 index f040258..81b6bc7 100644 --- a/src/QuAcK/print_qsGW.f90 +++ b/src/QuAcK/print_qsGW.f90 @@ -66,8 +66,8 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig write(*,'(2X,A27,F15.6)') 'qsGW GM total energy =',EqsGW + ENuc + EcGM write(*,'(2X,A27,F15.6)') 'qsGW exchange energy =',Ex write(*,'(2X,A27,F15.6)') 'qsGW correlation energy =',Ec - write(*,'(2X,A27,F15.6)') 'RPA correlation energy =',EcRPA - write(*,'(2X,A27,F15.6)') 'GM correlation energy =',EcGM + write(*,'(2X,A27,F15.6)') 'RPA@qsGW correlation energy =',EcRPA + write(*,'(2X,A27,F15.6)') 'GM@qsGW correlation energy =',EcGM write(*,*)'-------------------------------------------' write(*,*) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 06ee561..138cc0e 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,5 +1,5 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ENuc,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF) + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF) ! Compute linear response @@ -21,6 +21,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani logical,intent(in) :: triplet_manifold integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: PHF(nBas,nBas) @@ -38,7 +39,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani integer :: nBasSq integer :: ispin integer :: n_diis - double precision :: EcRPA + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) double precision :: EcGM double precision :: Conv double precision :: rcond @@ -128,7 +130,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(.not. GW0 .or. nSCF == 0) then call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) endif @@ -190,7 +192,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Print results call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM) + call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM) ! Increment @@ -235,8 +237,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(singlet_manifold) then ispin = 1 + EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif @@ -245,16 +249,27 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(triplet_manifold) then ispin = 2 + EcBSE(ispin) = 0d0 + call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'BSE@qsGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + endif end subroutine qsGW diff --git a/src/utils/Makefile b/src/utils/Makefile index 1ce9b74..09a8268 100644 --- a/src/utils/Makefile +++ b/src/utils/Makefile @@ -9,7 +9,7 @@ AR = ar ifeq ($(DEBUG),1) FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant else -FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O2 +FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O3 endif LIBS = ../../lib/*.a diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index d5d5fc2..c62f78f 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -15,7 +15,7 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) integer :: nEl(nspin) integer :: mu,nu,la,si double precision :: Ov,Kin,Nuc,ERI - double precision :: rs,R + double precision :: lambda ! Output variables @@ -25,14 +25,9 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) debug = .false. - open(unit=1,file='input/sph') - read(1,*) - read(1,*) rs + lambda = 1d0 - R = sqrt(dble(sum(nEl(:))))/2d0*rs -! R = 1d0 - - print*, 'Scaling integrals by ',R + print*, 'Scaling integrals by ',lambda open(unit=8 ,file='int/Ov.dat') open(unit=9 ,file='int/Kin.dat') @@ -53,7 +48,6 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) T = 0d0 do read(9,*,end=9) mu,nu,Kin - T(mu,nu) = Kin/R**2 enddo 9 close(unit=9) @@ -76,7 +70,7 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) do read(11,*,end=11) mu,nu,la,si,ERI - ERI = ERI/R + ERI = lambda*ERI ! <12|34> G(mu,nu,la,si) = ERI ! <32|14>