sph ready

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
+

input/basis

@@ -1,9 +1,51 @@
-1   3
+1   25
-S   3   1.00
-     38.3600000              0.0238090        
-      5.7700000              0.1548910        
-      1.2400000              0.4699870        
 S   1   1.00
-      0.2976000              1.0000000        
+      1.0000000              1.0000000        
-P   1   1.00
+S   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        

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

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

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

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
+

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

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
 

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
+  integer                       :: nSCF
-  double precision              :: eps,Conv
+  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
 

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
 

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

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
+  R = sqrt(dble(sum(nEl(:))))/2d0*rs
- rs = 1d0
+! 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>