From a8d26d6e36f48deb75406de701e0a0d0d5408c99 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 23 Mar 2020 17:26:12 +0100 Subject: [PATCH] numgrid test --- examples/molecule.H2 | 2 +- input/basis | 12 +- input/dft | 2 +- input/methods | 4 +- input/weight | 12 +- src/QuAcK/QuAcK.f90 | 6 +- src/QuAcK/excitation_density_Tmatrix.f90 | 10 ++ src/QuAcK/linear_response_C_pp.f90 | 2 +- src/QuAcK/linear_response_D_pp.f90 | 2 +- src/QuAcK/soG0T0.f90 | 4 - src/eDFT/AO_values_grid.f90 | 8 +- src/eDFT/eDFT.f90 | 12 +- src/eDFT/quadrature_grid.f90 | 10 +- src/eDFT/test_explicit.f90 | 156 ----------------------- 14 files changed, 55 insertions(+), 187 deletions(-) delete mode 100644 src/eDFT/test_explicit.f90 diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 8b4dc85..81c624a 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 2.3 + H 0. 0. 1.4 diff --git a/input/basis b/input/basis index c08616f..1e6d894 100644 --- a/input/basis +++ b/input/basis @@ -1,4 +1,4 @@ -1 6 +1 9 S 4 1 234.0000000 0.0025870 2 35.1600000 0.0195330 @@ -8,9 +8,17 @@ S 1 1 0.6669000 1.0000000 S 1 1 0.2089000 1.0000000 +S 1 + 1 0.0513800 1.0000000 P 1 1 3.0440000 1.0000000 P 1 1 0.7580000 1.0000000 +P 1 + 1 0.1993000 1.0000000 D 1 - 1 1.9650000 1.0000000 + 1 1.9650000 1.0000000 +D 1 + 1 0.4592000 1.0000000 + + diff --git a/input/dft b/input/dft index 638a677..c810a7c 100644 --- a/input/dft +++ b/input/dft @@ -13,7 +13,7 @@ # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 1 RVWN5 + 1 RMFL20 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) diff --git a/input/methods b/input/methods index a8010f9..e126029 100644 --- a/input/methods +++ b/input/methods @@ -3,7 +3,7 @@ # MP2 MP3 MP2-F12 T F F # CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD - F T F F F F F + F F F F F F F # CIS RPA RPAx ppRPA ADC F F F F F # G0F2 evGF2 G0F3 evGF3 @@ -11,6 +11,6 @@ # G0W0 evGW qsGW F F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F diff --git a/input/weight b/input/weight index c08616f..1e6d894 100644 --- a/input/weight +++ b/input/weight @@ -1,4 +1,4 @@ -1 6 +1 9 S 4 1 234.0000000 0.0025870 2 35.1600000 0.0195330 @@ -8,9 +8,17 @@ S 1 1 0.6669000 1.0000000 S 1 1 0.2089000 1.0000000 +S 1 + 1 0.0513800 1.0000000 P 1 1 3.0440000 1.0000000 P 1 1 0.7580000 1.0000000 +P 1 + 1 0.1993000 1.0000000 D 1 - 1 1.9650000 1.0000000 + 1 1.9650000 1.0000000 +D 1 + 1 0.4592000 1.0000000 + + diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index fd7eacf..0374372 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -669,9 +669,9 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - 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 soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(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 f541da9..c1de99e 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -222,6 +222,16 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end do +! do p=nC+1,nBas-nR +! do i=nC+1,nO +! do ab=1,nVV +! print*,p,i,ab,rho1(p,i,ab) +! end do +! end do +! end do + +! call matout(nVV,nOO,Y1) + end if end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 index ae7d399..9c93daa 100644 --- a/src/QuAcK/linear_response_C_pp.f90 +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -25,7 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) ! Define the chemical potential eF = e(nO) + e(nO+1) - eF = 0d0 +! eF = 0d0 ! Build C matrix for the singlet manifold diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 index 002a6bc..93167af 100644 --- a/src/QuAcK/linear_response_D_pp.f90 +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -25,7 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) ! Define the chemical potential eF = e(nO) + e(nO+1) - eF = 0d0 +! eF = 0d0 ! Build the D matrix for the singlet manifold diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index bf57a4f..8fc0f0b 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -92,10 +92,6 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) X1(:,:),Y1(:,:),rho1(:,:,:), & X2(:,:),Y2(:,:),rho2(:,:,:)) - - - rho2(:,:,:) = 0d0 - !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- diff --git a/src/eDFT/AO_values_grid.f90 b/src/eDFT/AO_values_grid.f90 index a5a72a6..67b672a 100644 --- a/src/eDFT/AO_values_grid.f90 +++ b/src/eDFT/AO_values_grid.f90 @@ -9,17 +9,17 @@ subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,E ! Input variables integer,intent(in) :: nBas,nShell - double precision,intent(in) :: CenterShell(maxShell,3) + double precision,intent(in) :: CenterShell(maxShell,ncart) integer,intent(in) :: TotAngMomShell(maxShell) integer,intent(in) :: KShell(maxShell) double precision,intent(in) :: DShell(maxShell,maxK) double precision,intent(in) :: ExpShell(maxShell,maxK) - double precision,intent(in) :: root(3,nGrid) + double precision,intent(in) :: root(ncart,nGrid) integer,intent(in) :: nGrid ! Local variables - integer :: atot,nShellFunction,a(3) + integer :: atot,nShellFunction,a(ncart) integer,allocatable :: ShellFunction(:,:) double precision :: rASq,xA,yA,zA,norm_coeff,prim @@ -28,7 +28,7 @@ subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,E ! Output variables double precision,intent(out) :: AO(nBas,nGrid) - double precision,intent(out) :: dAO(3,nBas,nGrid) + double precision,intent(out) :: dAO(ncart,nBas,nGrid) ! Initialization diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 333a77e..410faa4 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -56,15 +56,11 @@ program eDFT write(*,*) '******************************************' write(*,*) -! Test numgrid - -! call test_numgrid() - !------------------------------------------------------------------------ ! Read input information !------------------------------------------------------------------------ -! Read number of atoms, number of electrons of the system +! Read number of atoms, number of electroes of the system ! nC = number of core orbitals ! nO = number of occupied orbitals ! nV = number of virtual orbitals (see below) @@ -107,7 +103,7 @@ program eDFT allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas), & X(nBas,nBas),ERI(nBas,nBas,nBas,nBas),F(nBas,nBas)) -! Read integrals + ! Read integrals call cpu_time(start_int) @@ -133,6 +129,10 @@ program eDFT allocate(root(ncart,nGrid),weight(nGrid)) call quadrature_grid(nRad,nAng,nGrid,root,weight) +! Test numgrid + +! call test_numgrid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,nRad,nAng,nGrid,root,weight) + !------------------------------------------------------------------------ ! Calculate AO values at grid points !------------------------------------------------------------------------ diff --git a/src/eDFT/quadrature_grid.f90 b/src/eDFT/quadrature_grid.f90 index 420e80a..437b5f7 100644 --- a/src/eDFT/quadrature_grid.f90 +++ b/src/eDFT/quadrature_grid.f90 @@ -7,7 +7,9 @@ subroutine quadrature_grid(nRad,nAng,nGrid,root,weight) ! Input variables - integer,intent(in) :: nRad,nAng,nGrid + integer,intent(in) :: nRad + integer,intent(in) :: nAng + integer,intent(in) :: nGrid ! Local variables @@ -20,12 +22,12 @@ subroutine quadrature_grid(nRad,nAng,nGrid,root,weight) ! Output variables - double precision,intent(out) :: root(3,nGrid) + double precision,intent(out) :: root(ncart,nGrid) double precision,intent(out) :: weight(nGrid) ! Memory allocation - allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng)) + allocate(Radius(nRad),RadWeight(nRad),XYZ(ncart,nAng),XYZWeight(nAng)) ! Findthe radial grid @@ -49,7 +51,7 @@ subroutine quadrature_grid(nRad,nAng,nGrid,root,weight) write(*,50) write(*,20) do j = 1,nAng - write(*,60) j,(XYZ(k,j),k=1,3), XYZWeight(j) + write(*,60) j,(XYZ(k,j),k=1,ncart), XYZWeight(j) end do write(*,20) diff --git a/src/eDFT/test_explicit.f90 b/src/eDFT/test_explicit.f90 deleted file mode 100644 index 87cd777..0000000 --- a/src/eDFT/test_explicit.f90 +++ /dev/null @@ -1,156 +0,0 @@ -subroutine test_numgrid - - ! in this test we compute the grid for O - ! in the presence of two H in the H2O molecule - - use numgrid - use, intrinsic :: iso_c_binding, only: c_ptr - - implicit none - - real(8) :: radial_precision - integer :: min_num_angular_points - integer :: max_num_angular_points - integer :: num_centers - integer :: proton_charges(3) - real(8) :: x_coordinates_bohr(3) - real(8) :: y_coordinates_bohr(3) - real(8) :: z_coordinates_bohr(3) - real(8) :: alpha_max - integer :: max_l_quantum_number - real(8) :: alpha_min(3) - integer :: num_points - integer :: num_radial_points - integer :: num_angular_points - real(8), allocatable :: grid_x_bohr(:) - real(8), allocatable :: grid_y_bohr(:) - real(8), allocatable :: grid_z_bohr(:) - real(8), allocatable :: grid_w(:) - real(8), allocatable :: angular_grid_x_bohr(:) - real(8), allocatable :: angular_grid_y_bohr(:) - real(8), allocatable :: angular_grid_z_bohr(:) - real(8), allocatable :: angular_grid_w(:) - real(8), allocatable :: radial_grid_r_bohr(:) - real(8), allocatable :: radial_grid_w(:) - integer :: center_index - integer, parameter :: io_unit = 13 - real(8) :: ref(4), own(4) - integer :: ipoint - integer :: j - real(8) :: error - type(c_ptr) :: context - - radial_precision = 1.0d-12 - min_num_angular_points = 86 - max_num_angular_points = 302 - - num_centers = 3 - - proton_charges(1) = 8 - proton_charges(2) = 1 - proton_charges(3) = 1 - - x_coordinates_bohr(1) = 0.0d0 - x_coordinates_bohr(2) = 1.43d0 - x_coordinates_bohr(3) = -1.43d0 - - y_coordinates_bohr(1) = 0.0d0 - y_coordinates_bohr(2) = 0.0d0 - y_coordinates_bohr(3) = 0.0d0 - - z_coordinates_bohr(1) = 0.0d0 - z_coordinates_bohr(2) = 1.1d0 - z_coordinates_bohr(3) = 1.1d0 - - ! oxygen - alpha_max = 11720.0d0 - max_l_quantum_number = 2 - alpha_min(1) = 0.3023d0 - alpha_min(2) = 0.2753d0 - alpha_min(3) = 1.185d0 - - open(unit=io_unit, file='/Users/loos/Work/numgrid/test/reference_grid/cc-pVDZ.txt', access='sequential', action='read') - - context = numgrid_new_atom_grid(radial_precision, & - min_num_angular_points, & - max_num_angular_points, & - proton_charges(1), & - alpha_max, & - max_l_quantum_number, & - alpha_min) - - num_points = numgrid_get_num_grid_points(context) - if (num_points /= 16364) stop 1 - - allocate(grid_x_bohr(num_points)) - allocate(grid_y_bohr(num_points)) - allocate(grid_z_bohr(num_points)) - allocate(grid_w(num_points)) - - center_index = 0 - - call numgrid_get_grid(context, & - num_centers, & - center_index, & - x_coordinates_bohr, & - y_coordinates_bohr, & - z_coordinates_bohr, & - proton_charges, & - grid_x_bohr, & - grid_y_bohr, & - grid_z_bohr, & - grid_w) - - do ipoint = 1, num_points - read(io_unit, *) ref(1), ref(2), ref(3), ref(4) - own(1) = grid_x_bohr(ipoint) - own(2) = grid_y_bohr(ipoint) - own(3) = grid_z_bohr(ipoint) - own(4) = grid_w(ipoint) - do j = 1, 4 - error = own(j) - ref(j) - if (dabs(ref(j)) > 1.0e-15) error = error/ref(j) - if (dabs(error) > 1.0e-5) stop 1 - end do - end do - close(io_unit) - - deallocate(grid_x_bohr) - deallocate(grid_y_bohr) - deallocate(grid_z_bohr) - deallocate(grid_w) - - num_radial_points = numgrid_get_num_radial_grid_points(context) - if (num_radial_points /= 106) stop 1 - - allocate(radial_grid_r_bohr(num_radial_points)) - allocate(radial_grid_w(num_radial_points)) - - call numgrid_get_radial_grid(context, & - radial_grid_r_bohr, & - radial_grid_w) - - deallocate(radial_grid_r_bohr) - deallocate(radial_grid_w) - - num_angular_points = 14 - - allocate(angular_grid_x_bohr(num_angular_points)) - allocate(angular_grid_y_bohr(num_angular_points)) - allocate(angular_grid_z_bohr(num_angular_points)) - allocate(angular_grid_w(num_angular_points)) - - call numgrid_get_angular_grid(num_angular_points, & - angular_grid_x_bohr, & - angular_grid_y_bohr, & - angular_grid_z_bohr, & - angular_grid_w) - - deallocate(angular_grid_x_bohr) - deallocate(angular_grid_y_bohr) - deallocate(angular_grid_z_bohr) - deallocate(angular_grid_w) - - call numgrid_free_atom_grid(context) - -end subroutine