diff --git a/extract_sph.sh b/extract_sph.sh new file mode 100755 index 0000000..a833748 --- /dev/null +++ b/extract_sph.sh @@ -0,0 +1,86 @@ +#! /bin/bash + +for i in *.out +do + + echo + echo '***********************' + echo '*** ' $i ' ***' + 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 + echo '--- WFT information ---' + echo + + grep "Ec(MP2) =" $i + grep "Ec(CCSD) =" $i + grep "Ec(CCSD(T)) =" $i + + 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 + 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 + 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 + 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 + echo '*** DONE ***' + echo + +done + diff --git a/input/basis b/input/basis index b9ca7b5..e29e19f 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,51 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 25 S 1 1.00 - 0.2976000 1.0000000 -P 1 1.00 - 1.2750000 1.0000000 + 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 +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 diff --git a/input/methods b/input/methods index 187d53b..b2232ca 100644 --- a/input/methods +++ b/input/methods @@ -1,14 +1,14 @@ # RHF UHF MOM T F F # MP2 MP3 MP2-F12 - T F F + F F F # CCD CCSD CCSD(T) - T F F + F T T # CIS TDHF ADC - F F F + T T F # GF2 GF3 - F F + T F # G0W0 evGW qsGW - F F F + T T F # MCMP2 F diff --git a/input/molecule b/input/molecule index c78e87e..39a84f4 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 4 4 0 0 # Znuc x y z - He 0.0 0.0 0.0 + X 0.0 0.0 0.0 diff --git a/input/options b/input/options index 55951ff..e5dae4f 100644 --- a/input/options +++ b/input/options @@ -5,10 +5,10 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 F 1 # CIS/TDHF: singlet triplet - T F + T T # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 5 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize - 64 0.00001 T 5 F F F F F F F + 64 0.00001 T 5 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 new file mode 100755 index 0000000..f5af4fd --- /dev/null +++ b/run_sph.sh @@ -0,0 +1,39 @@ +#! /bin/bash + +Lmax=2 +Mmax=4 +rs=$1 + +if [ $# != 1 ] +then + echo "Please, specify rs value" +else + + echo "------------------------" + echo "Maxmium L value = " $Lmax + echo "Maxmium M value = " $Mmax + echo "------------------------" + echo + + for (( L=0 ; L<=$Lmax ; L++ )) ; do + + ne=$(bc -l <<< "(2*($L+1)*($L+1))") + echo + echo "------------------------" + echo "Number of electrons = " $ne + echo "------------------------" + echo + + for (( M=$L+1 ; M<=$Mmax ; M++ )) ; do + + 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}_${M}.out + + done + + done + +fi + diff --git a/src/QuAcK/CCSD.f90 b/src/QuAcK/CCSD.f90 index 01beb93..039d368 100644 --- a/src/QuAcK/CCSD.f90 +++ b/src/QuAcK/CCSD.f90 @@ -136,7 +136,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) call form_tau(nO,nV,t1,t2,tau) EcMP2 = 0.5d0*dot_product(pack(OOVV,.true.),pack(tau,.true.)) - write(*,'(1X,A10,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + write(*,'(1X,A20,1X,F10.6)') 'Ec(MP2) = ',EcMP2 ! Initialization @@ -223,6 +223,15 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) end if + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)' CCSD energy ' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A20,1X,F15.10)')' E(CCSD) = ',ECCSD + write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD) = ',EcCCSD + write(*,*)'----------------------------------------------------' + write(*,*) + ! Deallocate memory deallocate(hvv,hoo,hvo, & @@ -247,7 +256,7 @@ subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)' CCSDT(T) energy ' + write(*,*)' CCSD(T) energy ' write(*,*)'----------------------------------------------------' write(*,'(1X,A20,1X,F15.10)')' E(CCSD(T)) = ',ECCSD + EcCCT write(*,'(1X,A20,1X,F10.6)') ' Ec(CCSD(T)) = ',EcCCSD + EcCCT diff --git a/src/QuAcK/GF2.f90 b/src/QuAcK/GF2.f90 index 6eab83c..86e8c37 100644 --- a/src/QuAcK/GF2.f90 +++ b/src/QuAcK/GF2.f90 @@ -16,7 +16,9 @@ subroutine GF2(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0) ! Local variables integer :: nSCF,n_diis - double precision :: eps,Conv + double precision :: eps + double precision :: Conv + double precision :: rcond double precision,allocatable :: eGF2(:),eOld(:),Bpp(:,:,:),error_diis(:,:),e_diis(:,:) integer :: i,j,a,b,p,q @@ -99,7 +101,11 @@ subroutine GF2(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2) + +! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 eOld = eGF2 diff --git a/src/QuAcK/GF2_diag.f90 b/src/QuAcK/GF2_diag.f90 index 219a38c..7542fd4 100644 --- a/src/QuAcK/GF2_diag.f90 +++ b/src/QuAcK/GF2_diag.f90 @@ -15,8 +15,11 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0) ! Local variables - integer :: nSCF,n_diis - double precision :: eps,Conv + integer :: nSCF + integer :: n_diis + double precision :: eps + double precision :: Conv + double precision :: rcond double precision,allocatable :: eGF2(:),eOld(:),Bpp(:,:),error_diis(:,:),e_diis(:,:) integer :: i,j,a,b,p @@ -92,7 +95,9 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2) + + if(abs(rcond) < 1d-15) n_diis = 0 eOld = eGF2 diff --git a/src/QuAcK/GF3_diag.f90 b/src/QuAcK/GF3_diag.f90 index b6fbfd6..6f5b7cd 100644 --- a/src/QuAcK/GF3_diag.f90 +++ b/src/QuAcK/GF3_diag.f90 @@ -14,8 +14,10 @@ ! Local variables - integer :: nSCF,n_diis + integer :: nSCF + integer :: n_diis double precision :: eps,eps1,eps2,Conv + double precision :: rcond double precision,allocatable :: Sig2(:),SigInf(:),Sig3(:),eGF3(:),eOld(:) double precision,allocatable :: App(:,:),Bpp(:,:),Cpp(:,:),Dpp(:,:) double precision,allocatable :: Z(:),X2h1p(:),X1h2p(:),Sig2h1p(:),Sig1h2p(:) @@ -454,7 +456,9 @@ ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(nBas,nBas,n_diis,error_diis,e_diis,eGF3-eOld,eGF3) + call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF3-eOld,eGF3) + + if(abs(rcond) < 1d-15) n_diis = 0 ! Store result for next iteration diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 66ffa1f..134dfed 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -224,7 +224,9 @@ program QuAcK ! AO to MO integral transform for post-HF methods !------------------------------------------------------------------------ - call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis) +! call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis) + ERI_MO_basis = ERI_AO_basis + print*,'!!! MO = AO !!!' !------------------------------------------------------------------------ ! Compute MP2 energy diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index ceb175d..d5d5fc2 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 + double precision :: rs,R ! Output variables @@ -29,10 +29,10 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) read(1,*) read(1,*) rs -! rs = sqrt(dble(sum(nEl(:))))/2d0*rs - rs = 1d0 + R = sqrt(dble(sum(nEl(:))))/2d0*rs +! R = 1d0 - print*, 'Scaling integrals by ',rs + print*, 'Scaling integrals by ',R open(unit=8 ,file='int/Ov.dat') open(unit=9 ,file='int/Kin.dat') @@ -53,7 +53,7 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) T = 0d0 do read(9,*,end=9) mu,nu,Kin - T(mu,nu) = Kin/rs**2 + T(mu,nu) = Kin/R**2 enddo 9 close(unit=9) @@ -76,7 +76,7 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) do read(11,*,end=11) mu,nu,la,si,ERI - ERI = ERI/rs + ERI = ERI/R ! <12|34> G(mu,nu,la,si) = ERI ! <32|14>