numgrid/test/test.f90 at master · aakunitsa/numgrid · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
program test

   ! 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='../test/reference_grid.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 program