diff --git a/include/parameters.h b/include/parameters.h index 5aa6c3a..279fe1e 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -3,6 +3,7 @@ integer,parameter :: nsp = 3 integer,parameter :: maxEns = 10 integer,parameter :: maxShell = 512 + integer,parameter :: maxL = 7 integer,parameter :: n1eInt = 3 integer,parameter :: n2eInt = 4 integer,parameter :: n3eInt = 3 diff --git a/input/basis b/input/basis index 20bc6a7..fb05e68 100644 --- a/input/basis +++ b/input/basis @@ -1,49 +1,18 @@ -1 9 -S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 -S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 -S 1 - 1 0.0280500 1.0000000 -S 1 - 1 0.0086400 1.0000000 -P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 -P 1 - 1 0.0240300 1.0000000 -P 1 - 1 0.0057900 1.0000000 -D 1 - 1 0.1239000 1.0000000 -D 1 - 1 0.0725000 1.0000000 -2 5 +1 3 S 3 1 13.0100000 0.0196850 2 1.9620000 0.1379770 3 0.4446000 0.4781480 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 S 1 - 1 0.0297400 1.0000000 + 1 0.1220000 1.0000000 P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 1.0000000 - + 1 0.7270000 1.0000000 diff --git a/input/dft b/input/dft index c810a7c..86fa632 100644 --- a/input/dft +++ b/input/dft @@ -1,5 +1,5 @@ # Restricted or unrestricted KS calculation - LIM-RKS + GOK-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 diff --git a/input/methods b/input/methods index 93cf4d8..c1d2dda 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF MOM T F F # MP2 MP3 MP2-F12 - T F F + F F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD @@ -9,7 +9,7 @@ # CIS RPA RPAx ppRPA ADC F F F F F # G0F2 evGF2 G0F3 evGF3 - T F F F + F F F F # G0W0 evGW qsGW F F F # G0T0 evGT qsGT diff --git a/input/molecule b/input/molecule index e76247c..81c624a 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 2 2 0 0 + 2 1 1 0 0 # Znuc x y z - Li 0. 0. 0. - H 0. 0. 3.099 + H 0. 0. 0. + H 0. 0. 1.4 diff --git a/input/molecule.xyz b/input/molecule.xyz index fb94244..6edc99d 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - Li 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 1.6399202947 + H 0.0000000000 0.0000000000 0.0000000000 + H 0.0000000000 0.0000000000 0.7408481486 diff --git a/input/weight b/input/weight index 20bc6a7..fb05e68 100644 --- a/input/weight +++ b/input/weight @@ -1,49 +1,18 @@ -1 9 -S 8 - 1 1469.0000000 0.0007660 - 2 220.5000000 0.0058920 - 3 50.2600000 0.0296710 - 4 14.2400000 0.1091800 - 5 4.5810000 0.2827890 - 6 1.5800000 0.4531230 - 7 0.5640000 0.2747740 - 8 0.0734500 0.0097510 -S 8 - 1 1469.0000000 -0.0001200 - 2 220.5000000 -0.0009230 - 3 50.2600000 -0.0046890 - 4 14.2400000 -0.0176820 - 5 4.5810000 -0.0489020 - 6 1.5800000 -0.0960090 - 7 0.5640000 -0.1363800 - 8 0.0734500 0.5751020 -S 1 - 1 0.0280500 1.0000000 -S 1 - 1 0.0086400 1.0000000 -P 3 - 1 1.5340000 0.0227840 - 2 0.2749000 0.1391070 - 3 0.0736200 0.5003750 -P 1 - 1 0.0240300 1.0000000 -P 1 - 1 0.0057900 1.0000000 -D 1 - 1 0.1239000 1.0000000 -D 1 - 1 0.0725000 1.0000000 -2 5 +1 3 S 3 1 13.0100000 0.0196850 2 1.9620000 0.1379770 3 0.4446000 0.4781480 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 S 1 - 1 0.0297400 1.0000000 + 1 0.1220000 1.0000000 P 1 - 1 0.7270000 1.0000000 -P 1 - 1 0.1410000 1.0000000 - + 1 0.7270000 1.0000000 diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index b8b1dd7..f44ffb4 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -91,10 +91,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute linear response - call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR, & - nOOs,nVVs,eHF(:),ERI(:,:,:,:), & - Omega1s(:),X1s(:,:),Y1s(:,:), & - Omega2s(:),X2s(:,:),Y2s(:,:), & + call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:), & EcRPA(ispin)) EcRPA(ispin) = 1d0*EcRPA(ispin) @@ -105,8 +103,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute excitation densities for the T-matrix call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & - X1s(:,:),Y1s(:,:),rho1s(:,:,:), & - X2s(:,:),Y2s(:,:),rho2s(:,:,:)) + X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) !---------------------------------------------- ! Triplet manifold @@ -116,10 +113,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute linear response - call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR, & - nOOt,nVVt,eHF(:),ERI(:,:,:,:), & - Omega1t(:),X1t(:,:),Y1t(:,:), & - Omega2t(:),X2t(:,:),Y2t(:,:), & + call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:), & EcRPA(ispin)) EcRPA(ispin) = 3d0*EcRPA(ispin) @@ -130,16 +125,12 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute excitation densities for the T-matrix call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & - X1t(:,:),Y1t(:,:),rho1t(:,:,:), & - X2t(:,:),Y2t(:,:),rho2t(:,:,:)) + X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- - rho1s = 0d0 - rho1t = 0d0 - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 7a8349f..d8ad4dc 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -33,14 +33,21 @@ program QuAcK logical :: doXBS integer :: nShell - integer,allocatable :: TotAngMomShell(:),KShell(:) - double precision,allocatable :: CenterShell(:,:),DShell(:,:),ExpShell(:,:) + integer,allocatable :: TotAngMomShell(:) + integer,allocatable :: KShell(:) + double precision,allocatable :: CenterShell(:,:) + double precision,allocatable :: DShell(:,:) + double precision,allocatable :: ExpShell(:,:) + integer,allocatable :: max_ang_mom(:) + double precision,allocatable :: min_exponent(:,:) + double precision,allocatable :: max_exponent(:) integer :: TrialType double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:) double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:) - double precision,allocatable :: ERI_AO_basis(:,:,:,:),ERI_MO_basis(:,:,:,:) + double precision,allocatable :: ERI_AO_basis(:,:,:,:) + double precision,allocatable :: ERI_MO_basis(:,:,:,:) double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:) double precision :: start_QuAcK ,end_QuAcK ,t_QuAcK @@ -161,14 +168,16 @@ program QuAcK call read_geometry(nNuc,ZNuc,rNuc,ENuc) - allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), & - DShell(maxShell,maxK),ExpShell(maxShell,maxK)) + allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell),DShell(maxShell,maxK), & + ExpShell(maxShell,maxK),max_ang_mom(nNuc),min_exponent(nNuc,maxL+1),max_exponent(nNuc)) + !------------------------------------------------------------------------ ! Read basis set information !------------------------------------------------------------------------ - call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) + call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & + max_ang_mom,min_exponent,max_exponent) nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:)) !------------------------------------------------------------------------ @@ -671,9 +680,9 @@ program QuAcK call cpu_time(start_G0T0) call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF) - ! 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) + 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) call cpu_time(end_G0T0) diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index 5b176a7..1ac2d24 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -43,25 +43,17 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r cd = 0 do c=nO+1,nBas-nR - do d=c,nBas-nR + do d=c,nBas-nR cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + ERI(p,i,c,d)*X1(cd,ab) -! + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & -! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) - + rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,c,d)*X1(cd,ab) end do end do kl = 0 do k=nC+1,nO - do l=k,nO + do l=k,nO kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + ERI(p,i,k,l)*Y1(kl,ab) -! + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & -! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) - + rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,k,l)*Y1(kl,ab) end do end do @@ -73,30 +65,23 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r cd = 0 do c=nO+1,nBas-nR - do d=c,nBas-nR + do d=c,nBas-nR cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + ERI(p,nO+a,c,d)*X2(cd,ij) -! + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) & -! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d))) - + rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,c,d)*X2(cd,ij) end do end do kl = 0 do k=nC+1,nO - do l=k,nO + do l=k,nO kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + ERI(p,nO+a,k,l)*Y2(kl,ij) -! + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) & -! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l))) - + rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,k,l)*Y2(kl,ij) end do end do end do end do + end do end if @@ -114,20 +99,17 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r cd = 0 do c=nO+1,nBas-nR - do d=c+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) - + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) end do end do kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=k+1,nO kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) end do end do @@ -139,25 +121,23 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r cd = 0 do c=nO+1,nBas-nR - do d=c+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) - + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) end do end do kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=k+1,nO kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) end do end do end do end do + end do end if @@ -175,19 +155,17 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r cd = 0 do c=nO+1,nBas-nR - do d=c+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) end do end do kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=k+1,nO kl = kl + 1 - rho1(p,i,ab) = rho1(p,i,ab) & - + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) + rho1(p,i,ab) = rho1(p,i,ab) + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) end do end do @@ -201,17 +179,15 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do c=nO+1,nBas-nR do d=c+1,nBas-nR cd = cd + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij) end do end do kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=k+1,nO kl = kl + 1 - rho2(p,a,ij) = rho2(p,a,ij) & - + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) + rho2(p,a,ij) = rho2(p,a,ij) + (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij) end do end do diff --git a/src/QuAcK/quadrature_grid.f90 b/src/QuAcK/quadrature_grid.f90 index 420e80a..30a7152 100644 --- a/src/QuAcK/quadrature_grid.f90 +++ b/src/QuAcK/quadrature_grid.f90 @@ -27,7 +27,7 @@ subroutine quadrature_grid(nRad,nAng,nGrid,root,weight) allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng)) -! Findthe radial grid +! Find the radial grid scale = 1d0 call EulMac(Radius,RadWeight,nRad,scale) diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 index f1a3e00..2ad76c8 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -45,7 +45,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do p=nC+1,nBas-nR do i=nC+1,nO do cd=1,nVVs - eps = e(p) + e(i) - Omega1s(cd) + eps = e(p) + e(i) - Omega1s(cd) SigT(p) = SigT(p) + rho1s(p,i,cd)**2/eps enddo enddo @@ -56,7 +56,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do p=nC+1,nBas-nR do a=1,nV-nR do kl=1,nOOs - eps = e(p) + e(nO+a) - Omega2s(kl) + eps = e(p) + e(nO+a) - Omega2s(kl) SigT(p) = SigT(p) + rho2s(p,a,kl)**2/eps enddo enddo @@ -71,8 +71,8 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do p=nC+1,nBas-nR do i=nC+1,nO do cd=1,nVVt - eps = e(p) + e(i) - Omega1t(cd) - ! SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps + eps = e(p) + e(i) - Omega1t(cd) + SigT(p) = SigT(p) + rho1t(p,i,cd)**2/eps enddo enddo enddo @@ -82,8 +82,8 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do p=nC+1,nBas-nR do a=1,nV-nR do kl=1,nOOt - eps = e(p) + e(nO+a) - Omega2t(kl) - ! SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps + eps = e(p) + e(nO+a) - Omega2t(kl) + SigT(p) = SigT(p) + rho2t(p,a,kl)**2/eps enddo enddo enddo diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index 47e0998..fc218d4 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -79,8 +79,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) ! Compute linear response - call linear_response_pp(ispin,.true.,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV, & - seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2, & + call linear_response_pp(ispin,.true.,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:),sERI(:,:,:,:), & + Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:), & EcRPA) call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:)) @@ -109,8 +109,8 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) ! Solve the quasi-particle equation !---------------------------------------------- -! eG0T0(:) = seHF(:) + SigT(:) - eG0T0(:) = seHF(:) + Z(:)*SigT(:) + eG0T0(:) = seHF(:) + SigT(:) +! eG0T0(:) = seHF(:) + Z(:)*SigT(:) !---------------------------------------------- ! Dump results diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 410faa4..e627b46 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -17,6 +17,9 @@ program eDFT double precision,allocatable :: CenterShell(:,:) double precision,allocatable :: DShell(:,:) double precision,allocatable :: ExpShell(:,:) + integer,allocatable :: max_ang_mom(:) + double precision,allocatable :: min_exponent(:,:) + double precision,allocatable :: max_exponent(:) double precision,allocatable :: S(:,:) double precision,allocatable :: T(:,:) @@ -75,14 +78,15 @@ program eDFT call read_geometry(nNuc,ZNuc,rNuc,ENuc) - allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), & - DShell(maxShell,maxK),ExpShell(maxShell,maxK)) + allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell),DShell(maxShell,maxK), & + ExpShell(maxShell,maxK),max_ang_mom(nNuc),min_exponent(nNuc,maxL+1),max_exponent(nNuc)) !------------------------------------------------------------------------ ! Read basis set information !------------------------------------------------------------------------ - call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) + call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & + max_ang_mom,min_exponent,max_exponent) !------------------------------------------------------------------------ ! DFT options @@ -131,7 +135,7 @@ program eDFT ! Test numgrid -! call test_numgrid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,nRad,nAng,nGrid,root,weight) + call test_numgrid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent) !------------------------------------------------------------------------ ! Calculate AO values at grid points diff --git a/src/utils/read_basis.f90 b/src/utils/read_basis.f90 index 640b645..35add6e 100644 --- a/src/utils/read_basis.f90 +++ b/src/utils/read_basis.f90 @@ -1,4 +1,5 @@ -subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) +subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & + max_ang_mom,min_exponent,max_exponent) ! Read basis set information @@ -24,6 +25,10 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) integer,intent(out) :: nV + integer,intent(out) :: max_ang_mom(nNuc) + double precision,intent(out) :: min_exponent(nNuc,maxL+1) + double precision,intent(out) :: max_exponent(nNuc) + !------------------------------------------------------------------------ ! Primary basis set information !------------------------------------------------------------------------ @@ -37,57 +42,119 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh write(*,'(A28)') 'Gaussian basis set' write(*,'(A28)') '------------------' +! Initailization + nShell = 0 + max_ang_mom(:) = 0 + min_exponent(:,:) = huge(0d0) + max_exponent(:) = 0d0 + +!------------------------------------------------------------------------ +! Loop over atoms +!------------------------------------------------------------------------ do i=1,nNuc + read(2,*) iNuc,nShAt write(*,'(A28,1X,I16)') 'Atom n. ',iNuc write(*,'(A28,1X,I16)') 'number of shells ',nShAt write(*,'(A28)') '------------------' -! Basis function centers - +!------------------------------------------------------------------------ +! Loop over shells +!------------------------------------------------------------------------ do j=1,nShAt + nShell = nShell + 1 + + ! Basis function centers + do k=1,ncart CenterShell(nShell,k) = rNuc(iNuc,k) enddo -! Shell type and contraction degree + ! Shell type and contraction degree read(2,*) shelltype,KShell(nShell) - if(shelltype == "S") then + + select case (shelltype) + case ("S") + TotAngMomShell(nShell) = 0 write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) - elseif(shelltype == "P") then + + case ("P") + TotAngMomShell(nShell) = 1 - write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) - elseif(shelltype == "D") then + write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) + + case ("D") + TotAngMomShell(nShell) = 2 write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) - elseif(shelltype == "F") then + + case ("F") + TotAngMomShell(nShell) = 3 write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) - elseif(shelltype == "G") then + + case ("G") + TotAngMomShell(nShell) = 4 write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) - elseif(shelltype == "H") then + + case ("H") + TotAngMomShell(nShell) = 5 write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) - elseif(shelltype == "I") then + + case ("I") + TotAngMomShell(nShell) = 6 write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif + + case ("J") + + TotAngMomShell(nShell) = 7 + write(*,'(A28,1X,I16)') 'j-type shell with K = ',KShell(nShell) + + case default + + call print_warning('!!! Angular momentum too high !!!') + stop + + end select ! Read exponents and contraction coefficients - write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' - do k=1,Kshell(nShell) - read(2,*) kk,ExpShell(nShell,k),DShell(nShell,k) - write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo + write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' + do k=1,Kshell(nShell) + read(2,*) kk,ExpShell(nShell,k),DShell(nShell,k) + write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) + enddo + + min_exponent(iNuc,TotAngMomShell(nShell)+1) & + = min(min_exponent(iNuc,TotAngMomShell(nShell)+1),minval(ExpShell(nShell,1:KShell(nShell)))) + max_exponent(iNuc) = max(max_exponent(iNuc),maxval(ExpShell(nShell,:))) + max_ang_mom(iNuc) = max(max_ang_mom(iNuc),TotAngMomShell(nShell)) + enddo +!------------------------------------------------------------------------ +! End loop over shells +!------------------------------------------------------------------------ + write(*,'(A28)') '------------------' + +! print*,'maximum angular momemtum for atom n. ',iNuc,' = ' +! print*,max_ang_mom(iNuc) +! print*,'maximum exponent for atom n. ',iNuc,' = ' +! print*,max_exponent(iNuc) +! print*,'minimum exponent for atom n. ',iNuc,' = ' +! print*,min_exponent(iNuc,1:max_ang_mom(iNuc)+1) + enddo +!------------------------------------------------------------------------ +! End loop over atoms +!------------------------------------------------------------------------ ! Total number of shells