diff --git a/input/basis b/input/basis index fb05e68..6796e3b 100644 --- a/input/basis +++ b/input/basis @@ -1,18 +1,9 @@ 1 3 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 0.1220000 1.0000000 + 1 0.2976000 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.1220000 1.0000000 -P 1 - 1 0.7270000 1.0000000 + 1 1.2750000 1.0000000 diff --git a/input/molecule b/input/molecule index 81c624a..c78e87e 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 1 1 1 0 0 # Znuc x y z - H 0. 0. 0. - H 0. 0. 1.4 + He 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..797b5fc 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + He 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index fb05e68..6796e3b 100644 --- a/input/weight +++ b/input/weight @@ -1,18 +1,9 @@ 1 3 S 3 - 1 13.0100000 0.0196850 - 2 1.9620000 0.1379770 - 3 0.4446000 0.4781480 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 0.1220000 1.0000000 + 1 0.2976000 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.1220000 1.0000000 -P 1 - 1 0.7270000 1.0000000 + 1 1.2750000 1.0000000 diff --git a/src/eDFT/allocate_grid.f90 b/src/eDFT/allocate_grid.f90 new file mode 100644 index 0000000..f9e5627 --- /dev/null +++ b/src/eDFT/allocate_grid.f90 @@ -0,0 +1,67 @@ +subroutine allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent,nGrid) + +! Allocate quadrature grid with numgrid (Radovan Bast) + + use numgrid + use, intrinsic :: iso_c_binding, only: c_ptr + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + + integer,intent(in) :: nShell + integer,intent(in) :: TotAngMomShell(maxShell) + double precision,intent(in) :: ExpShell(maxShell,maxK) + + integer,intent(in) :: max_ang_mom(nNuc) + double precision,intent(in) :: min_exponent(nNuc,maxL+1) + double precision,intent(in) :: max_exponent(nNuc) + +! Local variables + + integer :: iNuc + + double precision :: radial_precision + integer :: min_num_angular_points + integer :: max_num_angular_points + integer :: num_points + + integer :: center_index + type(c_ptr) :: context + +! Output variables + + integer,intent(out) :: nGrid + +! Set useful variables + + radial_precision = 1d-12 + min_num_angular_points = 6 ! SG-0 + max_num_angular_points = 170 ! SG-3 + +! Get total number of grid points + + nGrid = 0 + + do iNuc=1,nNuc + + context = numgrid_new_atom_grid(radial_precision, & + min_num_angular_points, & + max_num_angular_points, & + int(ZNuc(iNuc)), & + max_exponent(iNuc), & + max_ang_mom(iNuc), & + min_exponent(iNuc,:) ) + + nGrid = nGrid + numgrid_get_num_grid_points(context) + + call numgrid_free_atom_grid(context) + + end do + +end subroutine allocate_grid diff --git a/src/eDFT/build_grid.f90 b/src/eDFT/build_grid.f90 new file mode 100644 index 0000000..26d5b00 --- /dev/null +++ b/src/eDFT/build_grid.f90 @@ -0,0 +1,104 @@ +subroutine build_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & + nGrid,weight,root) + +! Compute quadrature grid with numgrid (Radovan Bast) + + use numgrid + use, intrinsic :: iso_c_binding, only: c_ptr + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + + integer,intent(in) :: nShell + integer,intent(in) :: TotAngMomShell(maxShell) + double precision,intent(in) :: ExpShell(maxShell,maxK) + + integer,intent(in) :: max_ang_mom(nNuc) + double precision,intent(in) :: min_exponent(nNuc,maxL+1) + double precision,intent(in) :: max_exponent(nNuc) + + integer,intent(in) :: nGrid + +! Local variables + + integer :: iNuc + integer :: iG + + double precision :: radial_precision + integer :: min_num_angular_points + integer :: max_num_angular_points + integer :: num_points + + integer :: center_index + type(c_ptr) :: context + +! Output variables + + double precision,intent(out) :: root(ncart,nGrid) + double precision,intent(out) :: weight(nGrid) + +! Set useful variables + + radial_precision = 1d-12 + min_num_angular_points = 6 ! SG-0 + max_num_angular_points = 170 ! SG-3 + +!------------------------------------------------------------------------ +! Main loop over atoms +!------------------------------------------------------------------------ + + iG = 0 + + do iNuc=1,nNuc + + context = numgrid_new_atom_grid(radial_precision, & + min_num_angular_points, & + max_num_angular_points, & + int(ZNuc(iNuc)), & + max_exponent(iNuc), & + max_ang_mom(iNuc), & + min_exponent(iNuc,:) ) + + center_index = iNuc - 1 + num_points = numgrid_get_num_grid_points(context) + + call numgrid_get_grid(context, & + nNuc, & + center_index, & + rNuc(:,1), & + rNuc(:,2), & + rNuc(:,3), & + int(ZNuc(:)), & + root(1,iG+1:num_points), & + root(3,iG+1:num_points), & + root(3,iG+1:num_points), & + weight(iG+1:num_points) ) + + iG = iG + num_points + + call numgrid_free_atom_grid(context) + + end do + +!------------------------------------------------------------------------ +! End main loop over atoms +!------------------------------------------------------------------------ + +! Print grid + + write(*,*) ' ***********************' + write(*,*) ' *** QUADRATURE GRID ***' + write(*,*) ' ***********************' + write(*,'(A10,3X,3A15)') 'Grid point','X','Y','Z' + do iG=1,nGrid + write(*,'(I10,3X,4F15.10)') iG,weight(iG),root(:,iG) + end do + write(*,*) + +end subroutine build_grid diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index e627b46..9379d63 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -128,14 +128,16 @@ program eDFT !------------------------------------------------------------------------ ! Construct quadrature grid !------------------------------------------------------------------------ - call read_grid(SGn,nRad,nAng,nGrid) +! call read_grid(SGn,nRad,nAng,nGrid) + + call allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent,nGrid) allocate(root(ncart,nGrid),weight(nGrid)) - call quadrature_grid(nRad,nAng,nGrid,root,weight) -! Test numgrid + call build_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & + nGrid,weight,root) - call test_numgrid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent) +! call quadrature_grid(nRad,nAng,nGrid,root,weight) !------------------------------------------------------------------------ ! Calculate AO values at grid points