Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 2 : PROGRAM atom_upf_unittest
9 2 : USE atom_basis, ONLY: atom_basis_opt
10 : USE atom_kind_orbitals, ONLY: sgp_potential_conversion
11 : USE atom_operators, ONLY: atom_ppint_release,&
12 : atom_ppint_setup
13 : USE atom_sgp, ONLY: atom_sgp_potential_type,&
14 : atom_sgp_release,&
15 : fit_core_density_gaussians,&
16 : sgp_construction
17 : USE atom_types, ONLY: &
18 : CGTO_BASIS, atom_basis_gridrep, atom_basis_type, atom_integrals, atom_potential_type, &
19 : init_atom_basis, init_atom_basis_default_pp, init_atom_potential, release_atom_basis, &
20 : release_atom_potential, sgp_pseudo, upf_pseudo
21 : USE atom_upf, ONLY: atom_read_upf,&
22 : atom_release_upf,&
23 : atom_upfpot_type
24 : USE atom_utils, ONLY: atom_core_density,&
25 : atom_density
26 : USE cp_files, ONLY: close_file,&
27 : open_file
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: silent_print_level
32 : USE external_potential_types, ONLY: allocate_potential,&
33 : deallocate_potential,&
34 : init_potential,&
35 : sgp_potential_type
36 : USE ieee_arithmetic, ONLY: ieee_is_finite
37 : USE input_constants, ONLY: contracted_gto,&
38 : do_gapw_gct,&
39 : do_gapw_log,&
40 : do_rhf_atom
41 : USE input_cp2k_atom, ONLY: create_atom_section
42 : USE input_section_types, ONLY: section_release,&
43 : section_type,&
44 : section_vals_create,&
45 : section_vals_get_subs_vals,&
46 : section_vals_release,&
47 : section_vals_type,&
48 : section_vals_val_set
49 : USE kinds, ONLY: default_path_length,&
50 : dp
51 : USE libcp2k, ONLY: cp2k_finalize,&
52 : cp2k_init
53 : USE mathconstants, ONLY: fourpi,&
54 : pi
55 : USE orbital_pointers, ONLY: co,&
56 : init_orbital_pointers,&
57 : nco,&
58 : nso
59 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
60 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
61 : create_grid_atom,&
62 : deallocate_grid_atom,&
63 : grid_atom_type
64 : USE upf_types, ONLY: upf_has_complete_partial_waves
65 :
66 : IMPLICIT NONE
67 : LOGICAL :: source_root_argument
68 :
69 : CHARACTER(LEN=default_path_length) :: filename
70 : INTEGER :: i, io_unit
71 : REAL(KIND=dp) :: error(6)
72 210 : TYPE(atom_upfpot_type) :: pot
73 : TYPE(atom_sgp_potential_type) :: sgp
74 :
75 2 : CALL cp2k_init()
76 2 : io_unit = cp_logger_get_default_io_unit()
77 2 : CALL check_core_derivatives()
78 2 : CALL check_projector_addition_theorem()
79 2 : CALL check_atomic_potential()
80 2 : CALL check_upf_units()
81 2 : CALL check_nlcc_fit()
82 2 : CALL check_basis_origin()
83 2 : CALL check_basis_derivatives()
84 2 : CALL check_radial_grid_endpoints()
85 2 : CALL check_contracted_basis_reader()
86 2 : CALL check_basis_optimization_export()
87 2 : CALL check_uspp_import()
88 2 : CALL check_uspp_partial_waves()
89 :
90 : ! The same temporary is reused for successive KIND sections in Quickstep.
91 2 : ALLOCATE (pot%rho_nlcc(3))
92 8 : pot%rho_nlcc = 1.0_dp
93 2 : pot%core_correction = .TRUE.
94 10 : pot%econf = [2, 4, 0, 0]
95 2 : pot%pname = "previous kind"
96 2 : CALL atom_release_upf(pot)
97 2 : IF (ALLOCATED(pot%rho_nlcc)) ERROR STOP "UPF release leaked the core density"
98 2 : IF (pot%core_correction) ERROR STOP "UPF release retained core-correction metadata"
99 10 : IF (ANY(pot%econf /= 0)) ERROR STOP "UPF release retained the previous electron configuration"
100 2 : IF (pot%pname /= "") ERROR STOP "UPF release retained the previous potential name"
101 2 : CALL atom_release_upf(pot)
102 :
103 2 : ALLOCATE (sgp%a_nlcc(2), sgp%c_nlcc(2))
104 2 : sgp%has_nlcc = .TRUE.
105 2 : sgp%n_nlcc = 2
106 2 : sgp%is_nonlocal(2) = .TRUE.
107 2 : CALL atom_sgp_release(sgp)
108 2 : IF (ASSOCIATED(sgp%a_nlcc) .OR. ASSOCIATED(sgp%c_nlcc)) ERROR STOP "SGP release leaked NLCC arrays"
109 2 : IF (sgp%has_nlcc .OR. sgp%n_nlcc /= 0) ERROR STOP "SGP release retained NLCC metadata"
110 14 : IF (ANY(sgp%is_nonlocal)) ERROR STOP "SGP release retained nonlocal channels"
111 2 : CALL atom_sgp_release(sgp)
112 :
113 : ! Optional external UPF files exercise the real reader and Gaussian transformation.
114 : ! This also permits auditing a complete SSSP library without bundling it in CP2K.
115 4 : DO i = 1, COMMAND_ARGUMENT_COUNT()
116 2 : CALL GET_COMMAND_ARGUMENT(i, filename)
117 : ! do_regtest.py passes the CP2K source root before any optional UPF audit files.
118 2 : IF (i == 1) THEN
119 2 : INQUIRE (FILE=TRIM(filename)//"/tests/UNIT_TESTS", EXIST=source_root_argument)
120 2 : IF (source_root_argument) CYCLE
121 : END IF
122 0 : CALL atom_read_upf(pot, TRIM(filename))
123 0 : IF (ABS(SUM(pot%econf) - pot%zion) > 1.e-10_dp) ERROR STOP "UPF valence count mismatch"
124 0 : IF (.NOT. ALL(ieee_is_finite(pot%vlocal))) ERROR STOP "Non-finite UPF local potential"
125 0 : IF (ANY(pot%r(2:) <= pot%r(:SIZE(pot%r) - 1))) ERROR STOP "UPF mesh is not increasing"
126 0 : IF (ANY(pot%rab <= 0.0_dp)) ERROR STOP "Non-positive UPF mesh weights"
127 0 : IF (ALLOCATED(pot%dion)) THEN
128 0 : IF (.NOT. ALL(ieee_is_finite(pot%dion))) ERROR STOP "Non-finite UPF projector coefficients"
129 0 : IF (.NOT. ALL(ieee_is_finite(pot%beta))) ERROR STOP "Non-finite UPF projectors"
130 0 : IF (MAXVAL(ABS(pot%dion - TRANSPOSE(pot%dion))) > 1.e-10_dp) ERROR STOP "Non-Hermitian UPF D matrix"
131 : END IF
132 0 : CALL sgp_construction(sgp, upf_pot=pot, error=error)
133 0 : IF (.NOT. ALL(ieee_is_finite(error))) ERROR STOP "Non-finite UPF transformation errors"
134 0 : IF (pot%core_correction) THEN
135 0 : IF (.NOT. ASSOCIATED(sgp%a_nlcc)) ERROR STOP "Missing transformed UPF core density"
136 0 : IF (.NOT. ALL(ieee_is_finite(sgp%c_nlcc))) ERROR STOP "Non-finite transformed UPF core density"
137 : END IF
138 0 : IF (io_unit > 0) WRITE (io_unit, '(A,1X,A,1X,6ES16.7)') "UPF transformation:", TRIM(filename), error
139 0 : CALL atom_sgp_release(sgp)
140 4 : CALL atom_release_upf(pot)
141 : END DO
142 :
143 2 : CALL cp2k_finalize()
144 :
145 : CONTAINS
146 :
147 : ! **************************************************************************************************
148 : !> \brief ...
149 : ! **************************************************************************************************
150 2 : SUBROUTINE check_uspp_import()
151 : CHARACTER(LEN=*), PARAMETER :: fixture = "uspp-augmentation-test.upf"
152 :
153 : INTEGER :: unit
154 210 : TYPE(atom_upfpot_type) :: pot
155 : TYPE(cp_logger_type), POINTER :: logger
156 :
157 2 : logger => cp_get_default_logger()
158 2 : IF (logger%para_env%is_source()) THEN
159 1 : OPEN (NEWUNIT=unit, FILE=fixture, STATUS='REPLACE', ACTION='WRITE')
160 : WRITE (unit, '(A)') &
161 1 : '<UPF version="2.0.1">', &
162 1 : '<PP_HEADER element="Pt" pseudo_type="USPP" is_ultrasoft="true"', &
163 1 : ' z_valence="16" mesh_size="3" number_of_proj="2" l_max="1" />', &
164 1 : '<PP_MESH mesh="3">', &
165 1 : '<PP_R>', &
166 1 : '0.5 1.0 1.5', &
167 1 : '</PP_R>', &
168 1 : '<PP_RAB>', &
169 1 : '0.5 0.5 0.5', &
170 1 : '</PP_RAB>', &
171 1 : '</PP_MESH>', &
172 1 : '<PP_LOCAL type="real" size="3">', &
173 1 : '6 4 2', &
174 1 : '</PP_LOCAL>', &
175 1 : '<PP_NONLOCAL>', &
176 1 : '<PP_BETA.1 type="real" size="3" index="1" angular_momentum="0">', &
177 1 : '1 4 9', &
178 1 : '</PP_BETA.1>', &
179 1 : '<PP_BETA.2 type="real" size="3" index="2" angular_momentum="1">', &
180 1 : '2 -1 3', &
181 1 : '</PP_BETA.2>', &
182 1 : '<PP_DIJ>', &
183 1 : '2 0 0 8', &
184 1 : '</PP_DIJ>', &
185 1 : '<PP_AUGMENTATION q_with_l="T" nqf="0" nqlc="3">', &
186 1 : '<PP_Q>', &
187 1 : '2 0 0 3', &
188 1 : '</PP_Q>', &
189 1 : '<PP_QIJL.2.2.0 first_index="2" second_index="2" composite_index="3" angular_momentum="0">', &
190 1 : '-3 4 5', &
191 1 : '</PP_QIJL.2.2.0>', &
192 1 : '<PP_QIJL.1.1.0 type="real" size="3" first_index="1" second_index="1"', &
193 1 : ' composite_index="1" angular_momentum="0" is_null="false">', &
194 1 : '1 2 3', &
195 1 : '</PP_QIJL.1.1.0>', &
196 1 : '<PP_QIJL.1.2.1 first_index="1" second_index="2" composite_index="2" angular_momentum="1">', &
197 1 : '4 -5 6', &
198 1 : '</PP_QIJL.1.2.1>', &
199 1 : '<PP_QIJL.2.2.2 first_index="2" second_index="2" composite_index="3"', &
200 1 : ' angular_momentum="2" is_null="T" size="0"/>', &
201 1 : '</PP_AUGMENTATION>', &
202 1 : '</PP_NONLOCAL>', &
203 2 : '</UPF>'
204 1 : CLOSE (unit)
205 : END IF
206 2 : CALL logger%para_env%sync()
207 2 : CALL atom_read_upf(pot, fixture, read_augmentation=.TRUE.)
208 10 : IF (ANY(pot%econf(0:3) /= [1, 6, 9, 0])) ERROR STOP 'Wrong 62-electron core subtraction'
209 8 : IF (ANY(SHAPE(pot%qijl) /= [3, 3, 3]) .OR. LBOUND(pot%qijl, 3) /= 0) ERROR STOP 'Wrong augmentation shape'
210 14 : IF (MAXVAL(ABS(pot%qij - RESHAPE([2._dp, 0._dp, 0._dp, 3._dp], [2, 2]))) > 1.e-14_dp) THEN
211 0 : ERROR STOP 'PP_Q normalization changed'
212 : END IF
213 8 : IF (MAXVAL(ABS(pot%qijl(:, 1, 0) - [1._dp, 2._dp, 3._dp])) > 1.e-14_dp) ERROR STOP 'Wrong s-s augmentation'
214 8 : IF (MAXVAL(ABS(pot%qijl(:, 2, 1) - [4._dp, -5._dp, 6._dp])) > 1.e-14_dp) ERROR STOP 'Wrong s-p augmentation'
215 8 : IF (MAXVAL(ABS(pot%qijl(:, 3, 0) - [-3._dp, 4._dp, 5._dp])) > 1.e-14_dp) ERROR STOP 'Wrong p-p augmentation'
216 26 : IF (ANY(pot%qijl(:, :, 2) /= 0._dp)) ERROR STOP 'Null augmentation was not zero initialized'
217 8 : IF (MAXVAL(ABS(pot%beta(:, 2) - [4._dp, -1._dp, 2._dp])) > 1.e-14_dp) ERROR STOP 'USPP beta normalization changed'
218 2 : IF (ABS(pot%dion(2, 2) - 4._dp) > 1.e-14_dp) ERROR STOP 'USPP D is not in Hartree'
219 2 : CALL atom_release_upf(pot)
220 2 : IF (ALLOCATED(pot%qij) .OR. ALLOCATED(pot%qijl)) ERROR STOP 'Augmentation arrays were not released'
221 2 : CALL logger%para_env%sync()
222 2 : IF (logger%para_env%is_source()) THEN
223 1 : OPEN (NEWUNIT=unit, FILE=fixture, STATUS='OLD')
224 1 : CLOSE (unit, STATUS='DELETE')
225 : END IF
226 :
227 4 : END SUBROUTINE check_uspp_import
228 :
229 : ! **************************************************************************************************
230 : !> \brief ...
231 : ! **************************************************************************************************
232 2 : SUBROUTINE check_uspp_partial_waves()
233 : CHARACTER(LEN=*), PARAMETER :: fixture = 'uspp-partial-waves-test.upf'
234 :
235 : CHARACTER(LEN=5) :: has_wfc
236 : INTEGER :: unit, variant
237 210 : TYPE(atom_upfpot_type) :: potential
238 : TYPE(cp_logger_type), POINTER :: logger
239 :
240 2 : logger => cp_get_default_logger()
241 6 : DO variant = 1, 2
242 4 : has_wfc = 'false'
243 4 : IF (variant == 1) has_wfc = 'true'
244 4 : IF (logger%para_env%is_source()) THEN
245 2 : OPEN (NEWUNIT=unit, FILE=fixture, STATUS='REPLACE', ACTION='WRITE')
246 : WRITE (unit, '(A)') &
247 2 : '<UPF version="2.0.1">', &
248 2 : '<PP_HEADER element="H" pseudo_type="US" is_ultrasoft="true" is_paw="false"', &
249 2 : ' has_wfc="'//TRIM(has_wfc)//'" z_valence="1" mesh_size="3"', &
250 2 : ' number_of_proj="2" number_of_wfc="1" l_max="1"/>', &
251 2 : '<PP_MESH mesh="3">', &
252 2 : '<PP_R>', '0.5 1 1.5', '</PP_R>', &
253 2 : '<PP_RAB>', '0.5 0.5 0.5', '</PP_RAB>', '</PP_MESH>', &
254 2 : '<PP_LOCAL>', '6 4 2', '</PP_LOCAL>', &
255 2 : '<PP_NONLOCAL>', &
256 2 : '<PP_BETA.1 index="1" angular_momentum="0" size="3">', '1 4 9', '</PP_BETA.1>', &
257 2 : '<PP_BETA.2 index="2" angular_momentum="1" size="3">', '2 -1 3', '</PP_BETA.2>', &
258 2 : '<PP_DIJ>', '2 0 0 8', '</PP_DIJ>', &
259 2 : '<PP_AUGMENTATION q_with_l="T" nqf="0" nqlc="3">', &
260 2 : '<PP_Q>', '2 0 0 3', '</PP_Q>', &
261 2 : '<PP_QIJL.1.1.0 first_index="1" second_index="1" composite_index="1" angular_momentum="0">', &
262 2 : '1 2 3', '</PP_QIJL.1.1.0>', &
263 2 : '<PP_QIJL.1.2.1 first_index="1" second_index="2" composite_index="2" angular_momentum="1">', &
264 2 : '4 -5 6', '</PP_QIJL.1.2.1>', &
265 2 : '<PP_QIJL.2.2.0 first_index="2" second_index="2" composite_index="3" angular_momentum="0">', &
266 2 : '-3 4 5', '</PP_QIJL.2.2.0>', &
267 2 : '<PP_QIJL.2.2.2 first_index="2" second_index="2" composite_index="3" angular_momentum="2">', &
268 2 : '0 0 0', '</PP_QIJL.2.2.2>', '</PP_AUGMENTATION>', '</PP_NONLOCAL>', &
269 4 : '<PP_PSWFC>', '<PP_CHI.1 l="0">', '999 999 999', '</PP_CHI.1>', '</PP_PSWFC>'
270 2 : IF (variant == 1) THEN
271 : ! Reordered, deliberately unnormalized waves distinguish every index and side.
272 : WRITE (unit, '(A)') &
273 1 : '<PP_FULL_WFC number_of_wfc="2">', &
274 1 : '<PP_PSWFC.2 type="real" size="3" l="1">', '2 3 -4', '</PP_PSWFC.2>', &
275 1 : '<PP_AEWFC.1 index="1" l="0">', '-3 5 -7', '</PP_AEWFC.1>', &
276 1 : '<PP_AEWFC.2 index="2" l="1">', '5 -7 9', '</PP_AEWFC.2>', &
277 2 : '<PP_PSWFC.1 l="0">', '1 -2 3', '</PP_PSWFC.1>', '</PP_FULL_WFC>'
278 : END IF
279 2 : WRITE (unit, '(A)') '</UPF>'
280 2 : CLOSE (unit)
281 : END IF
282 4 : CALL logger%para_env%sync()
283 4 : IF (variant == 1) THEN
284 2 : CALL atom_read_upf(potential, fixture, read_header=.TRUE.)
285 2 : IF (ALLOCATED(potential%partial_waves)) ERROR STOP 'USPP full waves imported without augmentation request'
286 2 : CALL atom_release_upf(potential)
287 : END IF
288 4 : CALL atom_read_upf(potential, fixture, read_augmentation=.TRUE.)
289 4 : IF (potential%is_paw) ERROR STOP 'USPP full waves were classified as PAW'
290 4 : IF (variant == 1) THEN
291 2 : IF (.NOT. ALLOCATED(potential%partial_waves)) ERROR STOP 'Missing optional USPP full waves'
292 : ASSOCIATE (waves => potential%partial_waves)
293 6 : IF (ANY(SHAPE(waves%ae_wfc) /= [3, 2])) ERROR STOP 'USPP full-wave count differs from beta count'
294 18 : IF (ANY(waves%ae_wfc /= RESHAPE([-3._dp, 5._dp, -7._dp, 5._dp, -7._dp, 9._dp], [3, 2]))) THEN
295 0 : ERROR STOP 'USPP AE full-wave order or original normalization changed'
296 : END IF
297 18 : IF (ANY(waves%ps_wfc /= RESHAPE([1._dp, -2._dp, 3._dp, 2._dp, 3._dp, -4._dp], [3, 2]))) THEN
298 0 : ERROR STOP 'USPP PS full waves changed or were confused with atomic guesses'
299 : END IF
300 6 : IF (ANY(waves%angular_momenta /= [0, 1])) ERROR STOP 'Wrong USPP full-wave angular momenta'
301 : END ASSOCIATE
302 2 : IF (.NOT. upf_has_complete_partial_waves(potential)) ERROR STOP 'Complete original full waves not recognized'
303 2 : DEALLOCATE (potential%partial_waves%ps_wfc)
304 2 : IF (upf_has_complete_partial_waves(potential)) ERROR STOP 'Incomplete full-wave pair accepted'
305 : ELSE
306 2 : IF (upf_has_complete_partial_waves(potential)) ERROR STOP 'Atomic guess waves accepted as complete full waves'
307 2 : IF (ALLOCATED(potential%partial_waves)) ERROR STOP 'USPP without full waves allocated partial-wave data'
308 : END IF
309 4 : CALL atom_release_upf(potential)
310 4 : IF (ALLOCATED(potential%partial_waves)) ERROR STOP 'USPP full-wave data not released'
311 6 : CALL logger%para_env%sync()
312 : END DO
313 2 : IF (logger%para_env%is_source()) THEN
314 1 : OPEN (NEWUNIT=unit, FILE=fixture, STATUS='OLD')
315 1 : CLOSE (unit, STATUS='DELETE')
316 : END IF
317 2 : END SUBROUTINE check_uspp_partial_waves
318 :
319 : ! **************************************************************************************************
320 : !> \brief Check the repaired radial maps at the size that exposes endpoint cancellation.
321 : ! **************************************************************************************************
322 2 : SUBROUTINE check_radial_grid_endpoints()
323 : INTEGER, PARAMETER :: n = 6400, types(2) = [do_gapw_gct, do_gapw_log]
324 : REAL(dp), PARAMETER :: inner(2) = [4.83458313724687745e-18_dp, 6.97482911687147684e-18_dp], &
325 : outer(2) = [2.06843066219244770e17_dp, 57.5213142095328827_dp]
326 :
327 : INTEGER :: i
328 : REAL(dp) :: error
329 : TYPE(grid_atom_type), POINTER :: grid
330 :
331 2 : NULLIFY (grid)
332 6 : DO i = 1, SIZE(types)
333 4 : CALL allocate_grid_atom(grid)
334 4 : CALL create_grid_atom(grid, n, 1, 2, 0, types(i))
335 51204 : IF (.NOT. ALL(ieee_is_finite(grid%rad)) .OR. .NOT. ALL(ieee_is_finite(grid%wr))) THEN
336 0 : ERROR STOP 'Nonfinite radial grid'
337 : END IF
338 51208 : IF (ANY(grid%rad <= 0.0_dp) .OR. ANY(grid%wr <= 0.0_dp)) ERROR STOP 'Nonpositive radial grid'
339 25600 : IF (ANY(grid%rad(:n - 1) <= grid%rad(2:))) ERROR STOP 'Unordered radial points'
340 : ! Independent analytic moments include the r**2 radial measure.
341 25604 : error = ABS(SUM(grid%wr*EXP(-grid%rad)) - 2.0_dp)
342 25604 : error = MAX(error, ABS(SUM(grid%wr*EXP(-grid%rad)/grid%rad) - 1.0_dp))
343 25604 : error = MAX(error, ABS(SUM(grid%wr*EXP(-grid%rad2)) - SQRT(pi)/4.0_dp))
344 4 : IF (error > 1.e-10_dp) ERROR STOP 'Radial moment error'
345 : ! Independent high-precision endpoints reject clipping rounded cosines.
346 4 : IF (ABS(grid%rad(n)/inner(i) - 1.0_dp) > 1.e-13_dp .OR. &
347 0 : ABS(grid%rad(1)/outer(i) - 1.0_dp) > 1.e-13_dp) ERROR STOP 'Radial grid endpoint error'
348 6 : CALL deallocate_grid_atom(grid)
349 : END DO
350 2 : END SUBROUTINE check_radial_grid_endpoints
351 :
352 : ! **************************************************************************************************
353 : !> \brief Check primitive and contracted Gaussian derivatives on UPF meshes starting at zero.
354 : ! **************************************************************************************************
355 2 : SUBROUTINE check_basis_origin()
356 : INTEGER :: i, l, n, variant
357 : REAL(KIND=dp) :: expected(3), r(6001), rab(6001), scale
358 : TYPE(atom_basis_type) :: basis, radial
359 :
360 2 : CALL init_atom_basis_default_pp(basis)
361 14 : n = MAXVAL(basis%nbas)
362 : ! A long supplied mesh must not generate an unused singular auxiliary grid.
363 2 : r = [(0.002_dp*REAL(i - 1, dp), i=1, SIZE(r))]
364 12004 : rab = 0.002_dp
365 2 : CALL atom_basis_gridrep(basis, radial, r, rab)
366 24006 : IF (ANY(radial%grid%rad /= r) .OR. ANY(radial%grid%wr /= rab*r**2)) THEN
367 0 : ERROR STOP "Supplied UPF mesh changed during basis evaluation"
368 : END IF
369 12004 : IF (MAXVAL(ABS(radial%bf(:, 1, 0) - EXP(-basis%am(1, 0)*r**2))) > 1.e-14_dp) THEN
370 0 : ERROR STOP "Incorrect basis values on long UPF mesh"
371 : END IF
372 2 : CALL release_atom_basis(radial)
373 6 : DO variant = 1, 2
374 4 : scale = REAL(variant, dp)
375 4 : IF (variant == 2) THEN
376 2 : basis%basis_type = CGTO_BASIS
377 10 : ALLOCATE (basis%cm(n, n, 0:UBOUND(basis%am, 2)))
378 5054 : basis%cm = 0.0_dp
379 42 : DO i = 1, n
380 282 : basis%cm(i, i, :) = scale
381 : END DO
382 : END IF
383 4 : CALL atom_basis_gridrep(basis, radial, [0.0_dp, 1.e-16_dp, 0.1_dp], [0.1_dp, 0.1_dp, 0.1_dp])
384 5836 : IF (.NOT. ALL(ieee_is_finite(radial%bf)) .OR. .NOT. ALL(ieee_is_finite(radial%dbf)) .OR. &
385 0 : .NOT. ALL(ieee_is_finite(radial%ddbf))) ERROR STOP "Nonfinite Gaussian derivatives at the UPF origin"
386 32 : DO l = 0, UBOUND(basis%am, 2)
387 508 : DO i = 1, basis%nbas(l)
388 480 : expected = 0.0_dp
389 720 : IF (l == 0) expected = [1.0_dp, 0.0_dp, -2.0_dp*basis%am(i, l)]
390 480 : IF (l == 1) expected(2) = 1.0_dp
391 400 : IF (l == 2) expected(3) = 2.0_dp
392 1920 : IF (MAXVAL(ABS([radial%bf(1, i, l), radial%dbf(1, i, l), radial%ddbf(1, i, l)] - &
393 24 : scale*expected)) > 1.e-12_dp) ERROR STOP "Incorrect Gaussian origin limit"
394 : END DO
395 : END DO
396 6 : CALL release_atom_basis(radial)
397 : END DO
398 2 : CALL release_atom_basis(basis)
399 74 : END SUBROUTINE check_basis_origin
400 :
401 : ! **************************************************************************************************
402 : !> \brief Differentiate radial basis values and densities independently at finite radii.
403 : ! **************************************************************************************************
404 2 : SUBROUTINE check_basis_derivatives()
405 : REAL(dp), PARAMETER :: radius(4) = [0.08_dp, 0.3_dp, 0.9_dp, 2.1_dp], step = 2.e-4_dp
406 :
407 : INTEGER :: i, k, l, n, variant
408 : REAL(dp) :: derivative(5), errors(3), first, lap(5), &
409 : r(5), reference, rho(5), second
410 : REAL(dp), ALLOCATABLE :: pmat(:, :, :)
411 : TYPE(atom_basis_type) :: basis, radial
412 :
413 2 : CALL init_atom_basis_default_pp(basis)
414 14 : n = MAXVAL(basis%nbas)
415 : ! Nonunit exponents and signed mixing exercise primitives and genuine contractions.
416 14 : basis%nbas = 2
417 14 : basis%nprim = 2
418 14 : basis%am(1, :) = 0.37_dp
419 14 : basis%am(2, :) = 2.3_dp
420 10 : ALLOCATE (pmat(n, n, 0:UBOUND(basis%am, 2)))
421 2 : pmat = 0.0_dp
422 14 : pmat(1, 1, :) = 0.8_dp
423 14 : pmat(2, 2, :) = 0.4_dp
424 14 : pmat(1, 2, :) = -0.17_dp
425 14 : pmat(2, 1, :) = -0.17_dp
426 2 : errors = 0.0_dp
427 6 : DO variant = 1, 2
428 4 : IF (variant == 2) THEN
429 2 : basis%basis_type = CGTO_BASIS
430 10 : ALLOCATE (basis%cm(n, n, 0:UBOUND(basis%am, 2)))
431 5054 : basis%cm = 0.0_dp
432 14 : basis%cm(1, 1, :) = 0.7_dp
433 14 : basis%cm(2, 1, :) = -0.2_dp
434 14 : basis%cm(1, 2, :) = 0.1_dp
435 14 : basis%cm(2, 2, :) = 1.3_dp
436 : END IF
437 22 : DO k = 1, SIZE(radius)
438 96 : r = radius(k) + step*[-2.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 2.0_dp]
439 16 : CALL atom_basis_gridrep(basis, radial, r, SPREAD(step, 1, 5))
440 128 : DO l = 0, UBOUND(basis%am, 2)
441 304 : DO i = 1, 2
442 1152 : first = DOT_PRODUCT(radial%bf(:, i, l), [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp])/(12*step)
443 1152 : second = DOT_PRODUCT(radial%bf(:, i, l), [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp])/(12*step**2)
444 192 : errors(1) = MAX(errors(1), ABS(first - radial%dbf(3, i, l))/(1 + ABS(first)))
445 288 : errors(2) = MAX(errors(2), ABS(second - radial%ddbf(3, i, l))/(1 + ABS(second)))
446 : END DO
447 : END DO
448 16 : CALL atom_density(rho, pmat, radial, UBOUND(basis%am, 2), typ="RHO")
449 32 : CALL atom_density(derivative, pmat, radial, UBOUND(basis%am, 2), typ="DER")
450 112 : CALL atom_density(lap, pmat, radial, UBOUND(basis%am, 2), typ="LAP", rr=r*r)
451 96 : first = DOT_PRODUCT(rho, [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp])/(12*step)
452 96 : second = DOT_PRODUCT(rho, [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp])/(12*step**2)
453 16 : errors(1) = MAX(errors(1), ABS(first - derivative(3))/(1 + ABS(first)))
454 16 : reference = second + 2*first/radius(k)
455 16 : errors(3) = MAX(errors(3), ABS(reference - lap(3))/(1 + ABS(reference)))
456 20 : CALL release_atom_basis(radial)
457 : END DO
458 : END DO
459 2 : IF (io_unit > 0) WRITE (io_unit, '(A,3ES24.15)') "Atomic radial derivative/Laplacian scaled errors: ", errors
460 2 : DEALLOCATE (pmat)
461 2 : CALL release_atom_basis(basis)
462 8 : IF (MAXVAL(errors) > 2.e-7_dp) ERROR STOP "Atomic radial finite-difference derivatives disagree"
463 82 : END SUBROUTINE check_basis_derivatives
464 :
465 : ! **************************************************************************************************
466 : !> \brief Check fitted densities and their derivatives against an analytic reference.
467 : !> Evaluation points differ from the input mesh, and neither input mesh
468 : !> nor reference exponents coincide with the internal fitting sequence.
469 : ! **************************************************************************************************
470 2 : SUBROUTINE check_nlcc_fit()
471 : INTEGER, PARAMETER :: nr = 1601
472 : REAL(dp), PARAMETER :: amplitudes(3) = [0.02_dp, 1.7_dp, 0.4_dp], &
473 : exponents(3) = [0.15_dp, 2.3_dp, 35.0_dp]
474 :
475 : INTEGER :: i, mesh
476 : REAL(dp) :: density_error, derivative, derivative_error, dx, r(nr), rab(nr), radius, &
477 : reference, reference_derivative, rho(nr), value
478 2 : REAL(dp), ALLOCATABLE :: alpha(:), coeff(:)
479 :
480 6 : DO mesh = 1, 2
481 : dx = LOG(20.0_dp/1.e-6_dp)/REAL(nr - 1, dp)
482 6408 : DO i = 1, nr
483 6404 : IF (mesh == 1) THEN
484 3202 : r(i) = 20.0_dp*REAL(i - 1, dp)/REAL(nr - 1, dp)
485 3202 : rab(i) = 20.0_dp/REAL(nr - 1, dp)
486 : ELSE
487 3202 : r(i) = 1.e-6_dp*EXP(dx*REAL(i - 1, dp))
488 3202 : rab(i) = dx*r(i)
489 : END IF
490 25620 : rho(i) = SUM(amplitudes*EXP(-exponents*r(i)**2))
491 : END DO
492 4 : CALL fit_core_density_gaussians(r, rab, rho, alpha, coeff)
493 4 : density_error = 0.0_dp
494 4 : derivative_error = 0.0_dp
495 4808 : DO i = 0, 1200
496 4804 : radius = REAL(i, dp)*0.017_dp
497 19216 : reference = SUM(amplitudes*EXP(-exponents*radius**2))
498 19216 : reference_derivative = -2.0_dp*radius*SUM(exponents*amplitudes*EXP(-exponents*radius**2))
499 235396 : value = SUM(coeff*EXP(-alpha*radius**2))
500 235396 : derivative = -2.0_dp*radius*SUM(alpha*coeff*EXP(-alpha*radius**2))
501 4804 : density_error = MAX(density_error, ABS(value - reference))
502 4808 : derivative_error = MAX(derivative_error, ABS(derivative - reference_derivative))
503 : END DO
504 4 : IF (io_unit > 0) WRITE (io_unit, '(A,I2,2ES16.7)') &
505 2 : "Core Gaussian fit, density and derivative:", mesh, density_error, derivative_error
506 4 : IF (density_error > 2.e-7_dp) ERROR STOP "Gaussian core-density fit is inaccurate"
507 4 : IF (derivative_error > 3.e-6_dp) ERROR STOP "Gaussian core-density derivative is inaccurate"
508 4 : rho = 0.0_dp
509 4 : CALL fit_core_density_gaussians(r, rab, rho, alpha, coeff)
510 198 : IF (ANY(coeff /= 0.0_dp)) ERROR STOP "Zero core density generated nonzero Gaussian coefficients"
511 : END DO
512 2 : END SUBROUTINE check_nlcc_fit
513 :
514 : ! **************************************************************************************************
515 : !> \brief Read an independent miniature UPF fixture and check its numerical convention.
516 : !> The values are deliberately simple and are not a physical pseudopotential.
517 : !> Checking beta separately matters: an inverse rescaling of D leaves the NC
518 : !> operator unchanged but would require rescaling USPP augmentation data too.
519 : ! **************************************************************************************************
520 2 : SUBROUTINE check_upf_units()
521 : CHARACTER(LEN=*), PARAMETER :: fixture = &
522 : 'upf-units-test-with-a-potential-filename-longer-than-'// &
523 : 'the-default-eighty-character-string-buffer.upf'
524 :
525 : INTEGER :: unit
526 10758 : TYPE(atom_potential_type) :: atomic
527 : TYPE(cp_logger_type), POINTER :: logger
528 : TYPE(section_type), POINTER :: specification
529 : TYPE(section_vals_type), POINTER :: input, potential_section
530 :
531 2 : logger => cp_get_default_logger()
532 2 : IF (logger%para_env%is_source()) THEN
533 1 : CALL open_file(fixture, unit_number=unit, file_status="REPLACE", file_action="WRITE")
534 : WRITE (unit, '(A)') &
535 1 : '<UPF version="2.0.1">', &
536 1 : '<PP_HEADER element="He" pseudo_type="NC" z_valence="2"', &
537 1 : ' mesh_size="3" number_of_proj="2" l_max="0" />', &
538 1 : '<PP_MESH mesh="3">', &
539 1 : '<PP_R type="real" size="3">', &
540 1 : '0.5 1.0 1.5', &
541 1 : '</PP_R>', &
542 1 : '<PP_RAB type="real" size="3">', &
543 1 : '0.5 0.5 0.5', &
544 1 : '</PP_RAB>', &
545 1 : '</PP_MESH>', &
546 1 : '<PP_LOCAL type="real" size="3">', &
547 1 : '6.0 4.0 2.0', &
548 1 : '</PP_LOCAL>', &
549 1 : '<PP_NONLOCAL>', &
550 1 : '<PP_BETA.1 type="real" size="3" index="1" angular_momentum="0">', &
551 1 : '1.0 4.0 9.0', &
552 1 : '</PP_BETA.1>', &
553 1 : '<PP_BETA.2 type="real" size="3" index="2" angular_momentum="0">', &
554 1 : '2.0 -1.0 3.0', &
555 1 : '</PP_BETA.2>', &
556 1 : '<PP_DIJ type="real" size="4">', &
557 1 : '2.0 0.4 0.4 8.0', &
558 1 : '</PP_DIJ>', &
559 1 : '</PP_NONLOCAL>', &
560 2 : '</UPF>'
561 1 : CALL close_file(unit)
562 : END IF
563 2 : CALL logger%para_env%sync()
564 2 : NULLIFY (specification, input)
565 2 : CALL create_atom_section(specification)
566 2 : CALL section_vals_create(input, specification)
567 2 : potential_section => section_vals_get_subs_vals(input, 'POTENTIAL')
568 2 : CALL section_vals_val_set(potential_section, 'PSEUDO_TYPE', i_val=upf_pseudo)
569 2 : CALL section_vals_val_set(potential_section, 'POTENTIAL_FILE_NAME', c_val=fixture)
570 2 : CALL init_atom_potential(atomic, potential_section, 2)
571 : ASSOCIATE (potential => atomic%upf_pot)
572 8 : IF (MAXVAL(ABS(potential%vlocal - [3.0_dp, 2.0_dp, 1.0_dp])) > 1.e-14_dp) THEN
573 0 : ERROR STOP "UPF local potential is not in Hartree"
574 : END IF
575 8 : IF (MAXVAL(ABS(potential%beta(:, 1) - [2.0_dp, 4.0_dp, 6.0_dp])) > 1.e-14_dp) THEN
576 0 : ERROR STOP "UPF changed the first tabulated projector amplitude"
577 : END IF
578 8 : IF (MAXVAL(ABS(potential%beta(:, 2) - [4.0_dp, -1.0_dp, 2.0_dp])) > 1.e-14_dp) THEN
579 0 : ERROR STOP "UPF changed the second tabulated projector amplitude"
580 : END IF
581 16 : IF (MAXVAL(ABS(potential%dion - RESHAPE([1.0_dp, 0.2_dp, 0.2_dp, 4.0_dp], [2, 2]))) > 1.e-14_dp) THEN
582 0 : ERROR STOP "UPF D matrix is not in the tabulated-projector convention"
583 : END IF
584 : END ASSOCIATE
585 2 : CALL release_atom_potential(atomic)
586 2 : CALL section_vals_release(input)
587 2 : CALL section_release(specification)
588 2 : CALL logger%para_env%sync()
589 2 : IF (logger%para_env%is_source()) THEN
590 1 : CALL open_file(fixture, unit_number=unit, file_status="OLD")
591 1 : CALL close_file(unit, file_status="DELETE")
592 : END IF
593 2 : END SUBROUTINE check_upf_units
594 :
595 : ! **************************************************************************************************
596 : !> \brief Compare the converted atomic SGP operator with analytic Gaussian integrals.
597 : !> This checks that atomic guesses and the Quickstep Hamiltonian describe
598 : !> the same radial projectors and erf-regularized local Coulomb potential.
599 : ! **************************************************************************************************
600 2 : SUBROUTINE check_atomic_potential()
601 : INTEGER :: l
602 : REAL(KIND=dp) :: a, b, p, reference
603 : TYPE(atom_basis_type) :: basis
604 : TYPE(atom_integrals) :: integrals
605 10758 : TYPE(atom_potential_type) :: atomic
606 : TYPE(sgp_potential_type), POINTER :: potential
607 :
608 2 : NULLIFY (potential)
609 2 : CALL allocate_potential(potential)
610 2 : potential%zeff = 1.0_dp
611 2 : potential%core_charge_radius = 0.4_dp
612 2 : potential%alpha_core_charge = 0.5_dp/potential%core_charge_radius**2
613 2 : potential%n_nonlocal = 1
614 2 : potential%lmax = 2
615 8 : potential%is_nonlocal(0:2) = .TRUE.
616 2 : potential%n_local = 1
617 : ALLOCATE (potential%elec_conf(0:3), potential%a_nonlocal(1), &
618 : potential%c_nonlocal(1, 1, 0:2), potential%h_nonlocal(1, 0:2), &
619 2 : potential%a_local(1), potential%c_local(1))
620 10 : potential%elec_conf = [1, 0, 0, 0]
621 4 : potential%a_nonlocal = 0.7_dp
622 20 : potential%c_nonlocal = 1.0_dp/SQRT(fourpi)
623 14 : potential%h_nonlocal = 1.0_dp
624 4 : potential%a_local = 1.2_dp
625 4 : potential%c_local = 0.3_dp
626 2 : atomic%ppot_type = sgp_pseudo
627 2 : CALL sgp_potential_conversion(potential, atomic%sgp_pot)
628 2 : CALL init_atom_basis_default_pp(basis)
629 2 : CALL atom_ppint_setup(integrals, basis, atomic)
630 8 : DO l = 0, 2
631 6 : a = basis%am(8, l) + potential%a_nonlocal(1)
632 6 : reference = (GAMMA(l + 1.5_dp)/(2.0_dp*a**(l + 1.5_dp)))**2
633 8 : IF (ABS(integrals%hnl(8, 8, l) - reference) > 1.e-9_dp*ABS(reference)) THEN
634 0 : ERROR STOP "Atomic SGP nonlocal operator differs from analytic Gaussian integrals"
635 : END IF
636 : END DO
637 2 : p = 2.0_dp*basis%am(8, 0)
638 2 : b = 1.0_dp/(SQRT(2.0_dp)*potential%core_charge_radius)
639 : reference = -potential%zeff*b/(2.0_dp*p*SQRT(p + b*b)) + &
640 2 : potential%c_local(1)*GAMMA(1.5_dp)/(2.0_dp*(p + potential%a_local(1))**1.5_dp)
641 2 : IF (ABS(integrals%core(8, 8, 0) - reference) > 1.e-9_dp*ABS(reference)) THEN
642 0 : ERROR STOP "Atomic SGP local operator differs from analytic Gaussian integrals"
643 : END IF
644 2 : CALL atom_ppint_release(integrals)
645 2 : CALL release_atom_basis(basis)
646 2 : CALL deallocate_potential(potential)
647 464 : END SUBROUTINE check_atomic_potential
648 :
649 : ! **************************************************************************************************
650 : !> \brief Verify the physical angular kernel of Gaussian UPF projectors through l=5.
651 : !> For beta_l(r) = r**l exp(-a*r**2), stripping off the radial exponentials
652 : !> at unit vectors leaves sum_m Y_lm(u) Y_lm(v) = (2*l+1) P_l(u.v)/(4*pi).
653 : !> The reference uses the Legendre recurrence, independently of the Cartesian
654 : !> to spherical transformation used by the production projector integrals.
655 : ! **************************************************************************************************
656 2 : SUBROUTINE check_projector_addition_theorem()
657 : INTEGER, PARAMETER :: lmax = 5
658 :
659 : INTEGER :: i, j, l, lx, ly, lz, nc, ns
660 : REAL(KIND=dp) :: cross, diagonal, p(0:lmax), u(3), v(3), x
661 2 : REAL(KIND=dp), ALLOCATABLE :: pu(:), pv(:)
662 : TYPE(sgp_potential_type), POINTER :: potential
663 :
664 2 : NULLIFY (potential)
665 2 : CALL init_orbital_pointers(lmax)
666 2 : CALL init_spherical_harmonics(lmax, -1)
667 2 : CALL allocate_potential(potential)
668 2 : potential%n_nonlocal = 1
669 2 : potential%lmax = lmax
670 14 : potential%nppnl = SUM(nco(0:lmax))
671 2 : ALLOCATE (potential%c_nonlocal(1, 1, 0:lmax), potential%h_nonlocal(1, 0:lmax))
672 : ! atom_sgp stores the radial coefficients divided by sqrt(4*pi).
673 38 : potential%c_nonlocal = 1.0_dp/SQRT(fourpi)
674 26 : potential%h_nonlocal = 1.0_dp
675 2 : CALL init_potential(potential)
676 2 : u = [0.6_dp, 0.0_dp, 0.8_dp]
677 2 : v = [0.3_dp, -0.4_dp, 0.5_dp]
678 14 : v = v/NORM2(v)
679 8 : x = DOT_PRODUCT(u, v)
680 2 : p(0) = 1.0_dp
681 2 : p(1) = x
682 10 : DO l = 2, lmax
683 10 : p(l) = ((2*l - 1)*x*p(l - 1) - (l - 1)*p(l - 2))/l
684 : END DO
685 : nc = 0
686 : ns = 0
687 14 : DO l = 0, lmax
688 48 : ALLOCATE (pu(nso(l)), pv(nso(l)))
689 12 : pu = 0.0_dp
690 12 : pv = 0.0_dp
691 54 : DO lx = 0, l
692 166 : DO ly = 0, l - lx
693 112 : lz = l - lx - ly
694 112 : i = nc + co(lx, ly, lz)
695 1106 : DO j = 1, nso(l)
696 952 : pu(j) = pu(j) + potential%cprj_ppnl(i, ns + j)*u(1)**lx*u(2)**ly*u(3)**lz
697 1064 : pv(j) = pv(j) + potential%cprj_ppnl(i, ns + j)*v(1)**lx*v(2)**ly*v(3)**lz
698 : END DO
699 : END DO
700 : END DO
701 12 : diagonal = REAL(2*l + 1, dp)/fourpi
702 12 : cross = diagonal*p(l)
703 : IF (ABS(DOT_PRODUCT(pu, pu) - diagonal) > 1.e-12_dp .OR. &
704 228 : ABS(DOT_PRODUCT(pv, pv) - diagonal) > 1.e-12_dp .OR. &
705 : ABS(DOT_PRODUCT(pu, pv) - cross) > 1.e-12_dp) THEN
706 0 : ERROR STOP "SGP projector kernel violates the spherical-harmonic addition theorem"
707 : END IF
708 12 : nc = nc + nco(l)
709 12 : ns = ns + nso(l)
710 14 : DEALLOCATE (pu, pv)
711 : END DO
712 2 : CALL deallocate_potential(potential)
713 2 : END SUBROUTINE check_projector_addition_theorem
714 :
715 : ! **************************************************************************************************
716 : !> \brief Verify radial derivatives and the spherical Laplacian against finite differences.
717 : ! **************************************************************************************************
718 2 : SUBROUTINE check_core_derivatives()
719 : REAL(KIND=dp), PARAMETER :: h = 1.e-4_dp
720 :
721 : REAL(KIND=dp) :: deriv(4), lap(4), lap_fd(4), minus(4), &
722 : plus(4), r(4), rho(4)
723 10798 : TYPE(atom_potential_type) :: potential
724 :
725 2 : potential%ppot_type = sgp_pseudo
726 2 : potential%sgp_pot%has_nlcc = .TRUE.
727 2 : potential%sgp_pot%n_nlcc = 2
728 6 : potential%sgp_pot%a_nlcc(1:2) = [0.7_dp, 2.3_dp]
729 6 : potential%sgp_pot%c_nlcc(1:2) = [1.8_dp, -0.2_dp]
730 2 : r = [0.0_dp, 0.1_dp, 0.7_dp, 2.0_dp]
731 2 : rho = 0.0_dp
732 2 : plus = 0.0_dp
733 2 : minus = 0.0_dp
734 2 : deriv = 0.0_dp
735 2 : lap = 0.0_dp
736 2 : CALL atom_core_density(rho, potential, "RHO", r)
737 10 : CALL atom_core_density(plus, potential, "RHO", r + h)
738 10 : CALL atom_core_density(minus, potential, "RHO", r - h)
739 2 : CALL atom_core_density(deriv, potential, "DER", r)
740 2 : CALL atom_core_density(lap, potential, "LAP", r)
741 10 : IF (MAXVAL(ABS(deriv - (plus - minus)/(2.0_dp*h))) > 1.e-7_dp) THEN
742 0 : ERROR STOP "SGP core-density derivative differs from finite differences"
743 : END IF
744 10 : lap_fd = (plus - 2.0_dp*rho + minus)/h**2
745 2 : lap_fd(1) = 3.0_dp*lap_fd(1)
746 8 : lap_fd(2:) = lap_fd(2:) + (plus(2:) - minus(2:))/(h*r(2:))
747 10 : IF (MAXVAL(ABS(lap - lap_fd)) > 1.e-6_dp) ERROR STOP "SGP core-density Laplacian differs from finite differences"
748 2 : END SUBROUTINE check_core_derivatives
749 :
750 : ! **************************************************************************************************
751 : !> \brief Check independent shell, primitive and set counts in file and inline atomic bases.
752 : ! **************************************************************************************************
753 2 : SUBROUTINE check_contracted_basis_reader()
754 : INTEGER, PARAMETER :: contractions(3) = [11, 1, 1], lower_l(3) = [0, 1, 2], &
755 : primitives(3) = [11, 41, 2], sets(3) = [1, 1, 21], upper_l(3) = [3, 6, 3]
756 2 : CHARACTER(LEN=1600), ALLOCATABLE :: lines(:)
757 : INTEGER :: b, c, i, j, l, mode, n, nc, np, ns, row, &
758 : s, t, unit
759 : REAL(dp) :: alpha, coefficient, expected, &
760 : normalization
761 : TYPE(atom_basis_type) :: basis
762 : CHARACTER(LEN=*), PARAMETER :: fixture = "atom-contracted-basis-test"
763 :
764 : TYPE(cp_logger_type), POINTER :: logger
765 : TYPE(section_type), POINTER :: specification
766 : TYPE(section_vals_type), POINTER :: input, pp_basis
767 :
768 2 : logger => cp_get_default_logger()
769 8 : DO b = 1, SIZE(sets)
770 6 : np = primitives(b)
771 6 : nc = contractions(b)
772 6 : ns = sets(b)
773 18 : ALLOCATE (lines(1 + ns*(np + 1)))
774 6 : WRITE (lines(1), '(I0)') ns
775 6 : row = 1
776 52 : DO s = 1, ns
777 46 : row = row + 1
778 150 : WRITE (lines(row), '(*(I0,1X))') 1, lower_l(b), upper_l(b), np, (nc, l=lower_l(b), upper_l(b))
779 240 : DO i = 1, np
780 188 : row = row + 1
781 188 : alpha = REAL(10*s + i, dp)/20.0_dp
782 188 : WRITE (lines(row), '(*(F12.6,1X))') alpha, &
783 2610 : ((REAL(100*s + 10*i + c + 5*l, dp)/1000.0_dp, c=1, nc), l=lower_l(b), upper_l(b))
784 : END DO
785 : END DO
786 6 : IF (logger%para_env%is_source()) THEN
787 3 : OPEN (NEWUNIT=unit, FILE=fixture, STATUS="REPLACE", ACTION="WRITE")
788 123 : WRITE (unit, '(A)') 'H READER-TEST', (TRIM(lines(i)), i=1, SIZE(lines))
789 3 : CLOSE (unit)
790 : END IF
791 6 : CALL logger%para_env%sync()
792 18 : DO mode = 1, 2
793 12 : NULLIFY (specification, input)
794 12 : CALL create_atom_section(specification)
795 12 : CALL section_vals_create(input, specification)
796 12 : pp_basis => section_vals_get_subs_vals(input, 'PP_BASIS')
797 12 : CALL section_vals_val_set(pp_basis, 'BASIS_TYPE', i_val=contracted_gto)
798 12 : CALL section_vals_val_set(pp_basis, 'GRID_POINTS', i_val=32)
799 12 : CALL section_vals_val_set(pp_basis, 'BASIS_SET', c_val='READER-TEST')
800 12 : CALL section_vals_val_set(pp_basis, 'BASIS_SET_FILE_NAME', c_val=fixture)
801 12 : IF (mode == 2) THEN
802 246 : DO i = 1, SIZE(lines)
803 246 : CALL section_vals_val_set(pp_basis, 'BASIS%_DEFAULT_KEYWORD_', i_rep_val=i, c_val=TRIM(lines(i)))
804 : END DO
805 : END IF
806 12 : CALL init_atom_basis(basis, pp_basis, 1, 'PP')
807 84 : DO l = 0, UBOUND(basis%nbas, 1)
808 72 : n = 0
809 72 : IF (l >= lower_l(b) .AND. l <= upper_l(b)) n = ns
810 72 : IF (basis%nprim(l) /= n*np .OR. basis%nbas(l) /= n*nc) THEN
811 0 : ERROR STOP "Incorrect atomic contracted basis dimensions"
812 : END IF
813 288 : DO s = 1, n
814 1608 : DO i = 1, np
815 1332 : alpha = REAL(10*s + i, dp)/20.0_dp
816 1332 : j = (s - 1)*np + i
817 1332 : IF (ABS(basis%am(j, l) - alpha) > 1.e-14_dp) ERROR STOP "Incorrect contracted basis exponent"
818 : ! Inverse square root of the analytic radial Gaussian norm.
819 1332 : normalization = SQRT(2.0_dp*(2.0_dp*alpha)**(l + 1.5_dp)/GAMMA(l + 1.5_dp))
820 9588 : DO t = 1, n
821 19196 : DO c = 1, nc
822 9812 : coefficient = 0.0_dp
823 9812 : IF (s == t) coefficient = REAL(100*s + 10*i + c + 5*l, dp)/1000.0_dp
824 9812 : expected = normalization*coefficient
825 17864 : IF (ABS(basis%cm(j, (t - 1)*nc + c, l) - expected) > 1.e-12_dp*MAX(1.0_dp, ABS(expected))) THEN
826 0 : ERROR STOP "Incorrect contracted basis coefficient or normalization"
827 : END IF
828 : END DO
829 : END DO
830 : END DO
831 : END DO
832 : END DO
833 12 : CALL release_atom_basis(basis)
834 12 : CALL section_vals_release(input)
835 18 : CALL section_release(specification)
836 : END DO
837 6 : CALL logger%para_env%sync()
838 6 : IF (logger%para_env%is_source()) THEN
839 3 : OPEN (NEWUNIT=unit, FILE=fixture, STATUS="OLD")
840 3 : CLOSE (unit, STATUS="DELETE")
841 : END IF
842 8 : DEALLOCATE (lines)
843 : END DO
844 40 : END SUBROUTINE check_contracted_basis_reader
845 :
846 : ! **************************************************************************************************
847 : !> \brief Verify exported orbitals use the best basis restored after Powell optimization.
848 : ! **************************************************************************************************
849 2 : SUBROUTINE check_basis_optimization_export()
850 : CHARACTER(LEN=160) :: line
851 : INTEGER :: i, j, status, unit
852 : REAL(dp) :: alpha(2), coefficient(2), norm
853 2 : REAL(dp), POINTER :: exponents(:)
854 : TYPE(cp_logger_type), POINTER :: logger
855 : TYPE(section_type), POINTER :: specification
856 : TYPE(section_vals_type), POINTER :: input
857 :
858 2 : NULLIFY (specification, input)
859 2 : CALL create_atom_section(specification)
860 2 : CALL section_vals_create(input, specification)
861 2 : CALL section_vals_val_set(input, 'ELEMENT', c_val='He')
862 2 : CALL section_vals_val_set(input, 'ELECTRON_CONFIGURATION', c_val='1s2')
863 2 : CALL section_vals_val_set(input, 'MAX_ANGULAR_MOMENTUM', i_val=0)
864 2 : CALL section_vals_val_set(input, 'AE_BASIS%NUM_GTO', i_val=2)
865 2 : ALLOCATE (exponents(2))
866 6 : exponents = [2.0_dp, 0.4_dp]
867 2 : CALL section_vals_val_set(input, 'AE_BASIS%S_EXPONENTS', r_vals_ptr=exponents)
868 2 : CALL section_vals_val_set(input, 'METHOD%METHOD_TYPE', i_val=do_rhf_atom)
869 2 : CALL section_vals_val_set(input, 'OPTIMIZATION%EPS_SCF', r_val=1.e-10_dp)
870 2 : CALL section_vals_val_set(input, 'OPTIMIZATION%N_DIIS', i_val=8)
871 2 : CALL section_vals_val_set(input, 'OPTIMIZATION%MAX_ITER', i_val=100)
872 2 : CALL section_vals_val_set(input, 'POWELL%MAX_FUN', i_val=5)
873 2 : CALL section_vals_val_set(input, 'POWELL%STEP_SIZE', r_val=0.15_dp)
874 2 : CALL section_vals_val_set(input, 'POWELL%ACCURACY', r_val=1.e-5_dp)
875 2 : CALL section_vals_val_set(input, 'PRINT%FIT_BASIS%_SECTION_PARAMETERS_', i_val=silent_print_level)
876 2 : CALL atom_basis_opt(input)
877 2 : CALL section_vals_release(input)
878 2 : CALL section_release(specification)
879 :
880 2 : logger => cp_get_default_logger()
881 2 : CALL logger%para_env%sync()
882 2 : OPEN (NEWUNIT=unit, FILE='OPT_BASIS', STATUS='OLD', ACTION='READ')
883 : DO
884 12 : READ (unit, '(A)', IOSTAT=status) line
885 12 : IF (status /= 0) ERROR STOP 'Missing optimized basis orbital coefficients'
886 12 : IF (INDEX(line, 'L Quantum Number:') > 0) EXIT
887 : END DO
888 6 : DO i = 1, 2
889 6 : READ (unit, *) alpha(i), coefficient(i)
890 : END DO
891 2 : CLOSE (unit)
892 2 : IF (alpha(1) <= 2.0_dp) ERROR STOP 'Optimization did not restore the earlier best exponent'
893 : norm = 0.0_dp
894 6 : DO i = 1, 2
895 14 : DO j = 1, 2
896 : norm = norm + coefficient(i)*coefficient(j)* &
897 12 : (2.0_dp*SQRT(alpha(i)*alpha(j))/(alpha(i) + alpha(j)))**1.5_dp
898 : END DO
899 : END DO
900 : ! The export rounds its coefficients to eight decimal places.
901 2 : IF (.NOT. ieee_is_finite(norm) .OR. ABS(norm - 1.0_dp) > 1.e-7_dp) THEN
902 0 : ERROR STOP 'Exported orbital is inconsistent with the optimized basis'
903 : END IF
904 2 : CALL logger%para_env%sync()
905 2 : IF (logger%para_env%is_source()) THEN
906 1 : OPEN (NEWUNIT=unit, FILE='OPT_BASIS', STATUS='OLD')
907 1 : CLOSE (unit, STATUS='DELETE')
908 : END IF
909 2 : END SUBROUTINE check_basis_optimization_export
910 :
911 : END PROGRAM atom_upf_unittest
|
