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 upf_local_integrals_unittest
9 2 : USE ai_overlap_ppl, ONLY: ppl_integral,&
10 : ppl_integral_ri
11 : USE atom_upf, ONLY: atom_read_upf,&
12 : atom_release_upf,&
13 : atom_upfpot_type
14 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
15 : USE ieee_arithmetic, ONLY: ieee_is_finite
16 : USE kinds, ONLY: default_path_length,&
17 : dp
18 : USE libcp2k, ONLY: cp2k_finalize,&
19 : cp2k_init
20 : USE mathconstants, ONLY: pi
21 : USE orbital_pointers, ONLY: init_orbital_pointers,&
22 : ncoset
23 : USE upf_local_integrals, ONLY: upf_local_init,&
24 : upf_local_matrix,&
25 : upf_local_ri,&
26 : upf_local_type
27 :
28 : IMPLICIT NONE
29 : LOGICAL :: source_root_argument
30 :
31 : INTEGER, PARAMETER :: la = 4, lb = 4, na = 2, nb = 2
32 : REAL(KIND=dp), PARAMETER :: alpha(na) = [.1_dp, 3._dp], beta(nb) = [.3_dp, 2.1_dp], gamma(2) = [.2_dp, 4._dp]
33 : REAL(KIND=dp), PARAMETER :: coeff(1, 2) = RESHAPE([-.7_dp, 2.3_dp], [1, 2])
34 : INTEGER :: n, grid, geometry, i, j, k, ld, ma, mb, io_unit
35 : CHARACTER(len=default_path_length) :: filename
36 : REAL(KIND=dp) :: rabv(3), rac(3), rbc(3), dx, error, force_a(3), force_b(3), ref_a(3), ref_b(3), scale, threshold
37 4 : REAL(KIND=dp), ALLOCATABLE :: vri(:), dri(:, :), rri(:), rdri(:, :), energy_only(:, :)
38 4 : REAL(KIND=dp), ALLOCATABLE :: r(:), rab(:), potential(:), v(:, :), da(:, :, :), db(:, :, :), &
39 4 : reference(:, :), density(:, :), s(:, :, :), fs(:, :, :)
40 2 : TYPE(upf_local_type) :: cache
41 210 : TYPE(atom_upfpot_type) :: pot
42 2 : CALL cp2k_init()
43 2 : io_unit = cp_logger_get_default_io_unit()
44 2 : CALL init_orbital_pointers(2*MAX(la, lb) + 3)
45 10 : DO grid = 1, 4
46 8 : n = 2501
47 8 : IF (grid == 2) n = 2500
48 16 : ALLOCATE (r(n), rab(n), potential(n))
49 20014 : DO i = 1, n
50 20014 : IF (grid == 3) THEN
51 5002 : dx = LOG(1 + 12.5_dp/.001_dp)/(n - 1)
52 5002 : r(i) = .001_dp*(EXP((i - 1)*dx) - 1)
53 5002 : rab(i) = dx*(r(i) + .001_dp)
54 : ELSE
55 15004 : r(i) = (i - 1)*.005_dp
56 15004 : rab(i) = .005_dp
57 : END IF
58 : END DO
59 20014 : potential(:) = coeff(1, 1)*EXP(-gamma(1)*r*r) + coeff(1, 2)*EXP(-gamma(2)*r*r)
60 8 : ALLOCATE (pot%r(n), pot%rab(n), pot%vlocal(n))
61 20014 : pot%r(:) = r
62 20014 : pot%rab(:) = rab
63 8 : pot%zion = 6._dp
64 20014 : pot%vlocal(:) = potential
65 20014 : DO i = 1, n
66 20014 : IF (r(i) > 0._dp) THEN
67 19998 : pot%vlocal(i) = potential(i) - pot%zion*erf(SQRT(.8_dp)*r(i))/r(i)
68 : ELSE
69 8 : pot%vlocal(i) = potential(i) - 2*pot%zion*SQRT(.8_dp/pi)
70 : END IF
71 : END DO
72 8 : threshold = 0._dp
73 : ! The coarse default EPS_PPL may screen neighbors, never radial samples.
74 8 : IF (grid == 4) threshold = 1.e-2_dp
75 8 : CALL upf_local_init(pot, .8_dp, la + lb + 2, cache, threshold)
76 8 : ma = na*ncoset(la)
77 8 : mb = nb*ncoset(lb)
78 8 : ld = ncoset(2*MAX(la, lb) + 2)
79 176 : ALLOCATE (v(ma, mb), da(ma, mb, 3), db(ma, mb, 3), reference(ma, mb), density(ma, mb), s(ld, ld, la + lb + 3), fs(ld, ld, 4))
80 72 : ALLOCATE (vri(ma), dri(ma, 3), rri(ma), rdri(ma, 3), energy_only(ma, mb))
81 568 : DO j = 1, mb
82 39768 : DO i = 1, ma
83 39760 : density(i, j) = SIN(REAL(17*i + 3*j, dp))/REAL(ma*mb, dp)
84 : END DO
85 : END DO
86 40 : DO geometry = 1, 4
87 8 : SELECT CASE (geometry)
88 : CASE (1)
89 8 : rabv = 0; rac = 0
90 : CASE (2)
91 8 : rabv = [.3_dp, -.5_dp, .8_dp]; rac = [.2_dp, .4_dp, -.7_dp]
92 : CASE (3)
93 8 : rabv = [3._dp, 2._dp, -1._dp]; rac = [-1._dp, 2._dp, .3_dp]
94 : CASE (4)
95 32 : rabv = [-7._dp, 1._dp, 2._dp]; rac = [.5_dp, -.2_dp, 1.7_dp]
96 : END SELECT
97 128 : rbc = rac - rabv
98 32 : CALL upf_local_matrix(cache, la, alpha, lb, beta, rabv, rac, v, da, db)
99 32 : CALL upf_local_matrix(cache, la, alpha, lb, beta, rabv, rac, energy_only)
100 318144 : IF (MAXVAL(ABS(v - energy_only))/MAX(1._dp, MAXVAL(ABS(v))) > 1.e-12_dp) THEN
101 0 : ERROR STOP 'Energy-only and derivative local operators differ'
102 : END IF
103 32 : reference = 0; s = 0; fs = 0
104 : CALL ppl_integral(la, 0, na, SPREAD(100._dp, 1, na), alpha, lb, 0, nb, SPREAD(100._dp, 1, nb), beta, &
105 : 2, gamma, [1, 1], coeff, 100._dp, rabv, norm2(rabv), rac, norm2(rac), rbc, norm2(rbc), reference, s, &
106 320 : density, force_a, force_b, fs)
107 159072 : scale = MAX(1._dp, MAXVAL(ABS(reference)))
108 159072 : error = MAXVAL(ABS(v - reference))/scale
109 128 : DO k = 1, 3
110 477216 : ref_a(k) = SUM(density*da(:, :, k))
111 477248 : ref_b(k) = SUM(density*db(:, :, k))
112 : END DO
113 48 : IF (io_unit > 0) WRITE (io_unit, '(A,2I3,3ES18.9)') 'Native local matrix/derivative errors: ', grid, geometry, &
114 144 : error, MAXVAL(ABS(force_a - ref_a)), MAXVAL(ABS(force_b - ref_b))
115 32 : IF (error > 1.e-10_dp) ERROR STOP 'Local radial integral differs from analytic three-center reference'
116 224 : IF (MAXVAL(ABS(force_a - ref_a)) > 1.e-10_dp .OR. MAXVAL(ABS(force_b - ref_b)) > 1.e-10_dp) THEN
117 0 : ERROR STOP 'Local derivative mismatch'
118 : END IF
119 32 : CALL upf_local_ri(cache, la, alpha, rac, vri, dri)
120 32 : rri = 0; rdri = 0
121 : CALL ppl_integral_ri(la, 0, na, SPREAD(100._dp, 1, na), alpha, 2, gamma, [1, 1], coeff, 100._dp, &
122 224 : -rac, norm2(rac), rri, rdri)
123 9088 : scale = MAX(1._dp, MAXVAL(ABS(rri)), MAXVAL(ABS(rdri)))
124 9088 : error = MAX(MAXVAL(ABS(vri - rri)), MAXVAL(ABS(dri - rdri)))/scale
125 32 : IF (io_unit > 0) WRITE (io_unit, '(A,2I3,ES18.9)') 'Native local RI errors: ', grid, geometry, error
126 32 : IF (error > 1.e-10_dp) ERROR STOP 'Local RI operator differs from analytic Gaussian integrals'
127 32 : CALL upf_local_ri(cache, la, alpha, rac, vri)
128 2280 : IF (MAXVAL(ABS(vri - rri))/scale > 1.e-10_dp) ERROR STOP 'Energy-only RI local operator mismatch'
129 : END DO
130 8 : DEALLOCATE (r, rab, potential, v, da, db, reference, density, s, fs, vri, dri, rri, rdri, energy_only)
131 10 : CALL atom_release_upf(pot)
132 : END DO
133 : ! Optional real-library audit checks finiteness, not physical accuracy.
134 2 : ma = na*ncoset(la)
135 2 : mb = nb*ncoset(lb)
136 20 : ALLOCATE (v(ma, mb), da(ma, mb, 3), db(ma, mb, 3))
137 4 : DO i = 1, command_argument_count()
138 2 : CALL get_command_argument(i, filename)
139 : ! do_regtest.py passes the CP2K source root before any optional UPF audit files.
140 2 : IF (i == 1) THEN
141 2 : INQUIRE (FILE=TRIM(filename)//"/tests/UNIT_TESTS", EXIST=source_root_argument)
142 2 : IF (source_root_argument) CYCLE
143 : END IF
144 0 : CALL atom_read_upf(pot, TRIM(filename))
145 0 : CALL upf_local_init(pot, .8_dp, la + lb + 1, cache, 1.e-12_dp)
146 0 : rabv = [.3_dp, -.5_dp, .8_dp]; rac = [.2_dp, .4_dp, -.7_dp]
147 0 : CALL upf_local_matrix(cache, la, alpha, lb, beta, rabv, rac, v, da, db)
148 0 : IF (.NOT. ALL(ieee_is_finite(v)) .OR. .NOT. ALL(ieee_is_finite(da)) .OR. .NOT. ALL(ieee_is_finite(db))) THEN
149 0 : ERROR STOP 'Non-finite native local UPF integrals'
150 : END IF
151 0 : IF (io_unit > 0) WRITE (io_unit, '(A,1X,A,2ES18.9)') &
152 0 : 'Native local UPF:', TRIM(filename), cache%radial%radius, MAXVAL(ABS(v))
153 4 : CALL atom_release_upf(pot)
154 : END DO
155 2 : DEALLOCATE (v, da, db)
156 2 : CALL cp2k_finalize()
157 2 : END PROGRAM upf_local_integrals_unittest
|
