mirror of
https://github.com/QuantumPackage/qp2.git
synced 20241010 01:43:06 +02:00
72 lines
2.8 KiB
Fortran
72 lines
2.8 KiB
Fortran
subroutine example_becke_numerical_grid


implicit none


include 'constants.include.F'


BEGIN_DOC


! subroutine that illustrates the main features available in becke_numerical_grid


END_DOC


integer :: i,j,k,ipoint


double precision :: integral_1, integral_2,alpha,center(3)


print*,''


print*,'**************'


print*,'**************'


print*,'routine that illustrates the use of the grid'


print*,'**************'


print*,'This grid is built as the reunion of a spherical grid around each atom'


print*,'Each spherical grid contains a certain number of radial and angular points'


print*,''


print*,'n_points_integration_angular = ',n_points_integration_angular


print*,'n_points_radial_grid = ',n_points_radial_grid


print*,''


print*,'As an example of the use of the grid, we will compute the integral of a 3D gaussian'


! parameter of the gaussian: center of the gaussian is set to the first nucleus


center(1:3)=nucl_coord(1,1:3)


! alpha = exponent of the gaussian


alpha = 1.d0




print*,''


print*,'The first example uses the grid points as onedimensional array'


print*,'This is the mostly used representation of the grid'


print*,'It is the easyest way to use it with no drawback in terms of accuracy'


integral_1 = 0.d0


! you browse all the grid points as a onedimensional array


do i = 1, n_points_final_grid


double precision :: weight, r(3)


! you get x, y and z of the ith grid point


r(1) = final_grid_points(1,i)


r(2) = final_grid_points(2,i)


r(3) = final_grid_points(3,i)


weight = final_weight_at_r_vector(i)


double precision :: distance, f_r


! you compute the function to be integrated


distance = dsqrt( (r(1)  center(1))**2 + (r(2)  center(2))**2 + (r(3)  center(3))**2 )


f_r = dexp(alpha * distance*distance)


! you add the contribution of the grid point to the integral


integral_1 += f_r * weight


enddo


print*,'integral_1 =',integral_1


print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5


print*,''


print*,''


print*,'The second example uses the grid points as a collection of spherical grids centered on each atom'


print*,'This is mostly useful if one needs to split contributions between radial/angular/atomic of an integral'


! you browse the nuclei


do i = 1, nucl_num


! you browse the radial points attached to each nucleus


do j = 1, n_points_radial_grid


! you browse the angular points attached to each radial point of each nucleus


do k = 1, n_points_integration_angular


r(1) = grid_points_per_atom(1,k,j,i)


r(2) = grid_points_per_atom(2,k,j,i)


r(3) = grid_points_per_atom(3,k,j,i)


weight = final_weight_at_r(k,j,i)


distance = dsqrt( (r(1)  center(1))**2 + (r(2)  center(2))**2 + (r(3)  center(3))**2 )


f_r = dexp(alpha * distance*distance)


integral_2 += f_r * weight


enddo


enddo


enddo


print*,'integral_2 =',integral_2


print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5


print*,''


end
