diff --git a/GoSph b/GoSph index 8055886..282b774 100755 --- a/GoSph +++ b/GoSph @@ -9,9 +9,9 @@ if [ $# = 2 ] then cp examples/molecule.Sph_"$1" input/molecule cp examples/basis.Sph.Ylm"$2" input/basis -cp int/Sph_ERI_"$2".dat int/ERI.dat -cp int/Sph_Kin_"$2".dat int/Kin.dat -cp int/Sph_Nuc_"$2".dat int/Nuc.dat -cp int/Sph_Ov_"$2".dat int/Ov.dat +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 fi diff --git a/include/parameters.h b/include/parameters.h index 2109908..20d7562 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -2,7 +2,7 @@ integer,parameter :: nspin = 2 integer,parameter :: nsp = 3 integer,parameter :: maxEns = 10 - integer,parameter :: maxShell = 50 + integer,parameter :: maxShell = 1024 integer,parameter :: n1eInt = 3 integer,parameter :: n2eInt = 4 integer,parameter :: n3eInt = 3 diff --git a/input/basis b/input/basis index b9ca7b5..00205a4 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,129 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 64 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 +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 +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 40fd2a7..341e057 100644 --- a/input/methods +++ b/input/methods @@ -9,6 +9,6 @@ # GF2 GF3 F F # G0W0 evGW qsGW - T T F + T F F # MCMP2 F diff --git a/input/molecule b/input/molecule index c78e87e..8e32e16 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 25 25 0 0 # Znuc x y z - He 0.0 0.0 0.0 + X 0.0 0.0 0.0 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 7901fbf..66ffa1f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -165,7 +165,7 @@ program QuAcK ! Read integrals - call read_integrals(nBas,S,T,V,Hc,ERI_AO_basis) + call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis) ! Compute orthogonalization matrix diff --git a/src/QuAcK/RHF.f90 b/src/QuAcK/RHF.f90 index 387de91..871f800 100644 --- a/src/QuAcK/RHF.f90 +++ b/src/QuAcK/RHF.f90 @@ -16,12 +16,25 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH ! Local variables - integer :: nSCF,nBasSq,n_diis - double precision :: ET,EV,EJ,EK,Conv,Gap + integer :: nSCF + integer :: nBasSq + integer :: n_diis + double precision :: ET,EV,EJ,EK + double precision :: Conv + double precision :: Gap double precision :: rcond double precision,external :: trace_matrix - double precision,allocatable :: error(:,:),error_diis(:,:),F_diis(:,:) - double precision,allocatable :: J(:,:),K(:,:),cp(:,:),F(:,:),Fp(:,:) + double precision,allocatable :: error(:,:) + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: F(:,:) + double precision,allocatable :: Fp(:,:) + double precision,allocatable :: ON(:) + + integer :: i ! Output variables @@ -44,8 +57,8 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH ! Memory allocation - allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & - cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas), & + allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & + cp(nBas,nBas),Fp(nBas,nBas),F(nBas,nBas),ON(nBas), & error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Guess coefficients and eigenvalues @@ -65,6 +78,13 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) +! ON(:) = 0d0 +! do i=1,nO +! ON(i) = 1d0 +! ON(i) = dble(2*i-1) +! end do + +! call density_matrix(nBas,ON,c,P) ! Initialization @@ -123,6 +143,8 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERH P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) +! call density_matrix(nBas,ON,c,P) + ! Compute HF energy ERHF = trace_matrix(nBas,matmul(P,Hc)) & diff --git a/src/utils/read_integrals.f90 b/src/utils/read_integrals.f90 index e6f9dd6..3ea2ac7 100644 --- a/src/utils/read_integrals.f90 +++ b/src/utils/read_integrals.f90 @@ -1,8 +1,9 @@ -subroutine read_integrals(nBas,S,T,V,Hc,G) +subroutine read_integrals(nEl,nBas,S,T,V,Hc,G) ! Read one- and two-electron integrals from files implicit none + include 'parameters.h' ! Input variables @@ -11,9 +12,10 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) ! Local variables logical :: debug + integer :: nEl(nspin) integer :: mu,nu,la,si double precision :: Ov,Kin,Nuc,ERI - double precision :: scale + double precision :: rs ! Output variables @@ -23,7 +25,13 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) debug = .false. - scale = 1d0 + open(unit=1,file='input/sph') + read(1,*) + read(1,*) rs + + rs = sqrt(dble(sum(nEl(:))))/2d0*rs + + print*, 'Scaling integrals by ',rs open(unit=8 ,file='int/Ov.dat') open(unit=9 ,file='int/Kin.dat') @@ -44,7 +52,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) T = 0d0 do read(9,*,end=9) mu,nu,Kin - T(mu,nu) = Kin/scale**2 + T(mu,nu) = Kin/rs**2 enddo 9 close(unit=9) @@ -67,7 +75,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G) do read(11,*,end=11) mu,nu,la,si,ERI - ERI = ERI/scale + ERI = ERI/rs ! <12|34> G(mu,nu,la,si) = ERI ! <32|14>