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_projector_integrals_unittest
9 2 : USE ai_moments, ONLY: moment
10 : USE ai_overlap, ONLY: overlap
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 orbital_pointers, ONLY: nco,&
21 : ncoset
22 : USE orbital_transformation_matrices, ONLY: orbtramat
23 : USE upf_projector_integrals, ONLY: upf_projector_init,&
24 : upf_projector_overlap,&
25 : upf_projector_type
26 :
27 : IMPLICIT NONE
28 : LOGICAL :: source_root_argument
29 :
30 : INTEGER, PARAMETER :: la = 4, lb = 3, nder = 2, na = 5, nb = 5
31 : REAL(KIND=dp), PARAMETER :: alpha(na) = [0.04_dp, 0.3_dp, 1.2_dp, 5.0_dp, 24.0_dp]
32 : REAL(KIND=dp), PARAMETER :: eta(nb) = [0.3_dp, 0.7_dp, 1.1_dp, 2.2_dp, 0.4_dp]
33 : CHARACTER(LEN=default_path_length) :: filename
34 : INTEGER :: grid, nr, i, j, p, m, l, geometry, id, ld, ldai, first, shell, q, io_unit
35 : REAL(KIND=dp) :: dx, rac(3), origin(3), error, scale, expected_d
36 4 : REAL(KIND=dp), ALLOCATABLE :: sab(:, :), work(:, :, :), moments(:, :, :), integrals(:, :, :), reference(:, :, :)
37 210 : TYPE(atom_upfpot_type) :: pot
38 2 : TYPE(upf_projector_type) :: projectors
39 :
40 2 : CALL cp2k_init()
41 2 : io_unit = cp_logger_get_default_io_unit()
42 2 : pot%number_of_proj = nb
43 2 : ALLOCATE (pot%lbeta(nb))
44 12 : pot%lbeta(:) = [0, 1, 2, 3, 1]
45 2 : ALLOCATE (pot%dion(nb, nb))
46 62 : pot%dion = 0.0_dp
47 12 : DO i = 1, nb
48 12 : pot%dion(i, i) = 0.3_dp*i - 0.8_dp
49 : END DO
50 2 : pot%dion(2, 5) = 0.17_dp
51 2 : pot%dion(5, 2) = 0.17_dp
52 :
53 8 : DO grid = 1, 3
54 6 : nr = 2501
55 6 : IF (grid == 2) nr = 2500
56 12 : ALLOCATE (pot%r(nr), pot%rab(nr), pot%beta(nr, nb))
57 15010 : DO i = 1, nr
58 15010 : IF (grid == 3) THEN
59 5002 : dx = LOG(1.0_dp + 12.5_dp/0.001_dp)/(nr - 1)
60 5002 : pot%r(i) = 0.001_dp*(EXP((i - 1)*dx) - 1.0_dp)
61 5002 : pot%rab(i) = dx*(pot%r(i) + 0.001_dp)
62 : ELSE
63 10002 : pot%r(i) = (i - 1)*0.005_dp
64 10002 : pot%rab(i) = 0.005_dp
65 : END IF
66 : END DO
67 36 : DO i = 1, nb
68 75056 : pot%beta(:, i) = pot%r**pot%lbeta(i)*EXP(-eta(i)*pot%r**2)
69 : END DO
70 6 : CALL upf_projector_init(pot, la + nder, projectors)
71 : ! Independent expansion of D, including non-diagonal coupling of two p projectors.
72 6 : p = 0
73 36 : DO i = 1, nb
74 150 : DO m = -pot%lbeta(i), pot%lbeta(i)
75 114 : p = p + 1
76 114 : q = 0
77 714 : DO j = 1, nb
78 2850 : DO l = -pot%lbeta(j), pot%lbeta(j)
79 2166 : q = q + 1
80 2166 : expected_d = 0.0_dp
81 2166 : IF (pot%lbeta(i) == pot%lbeta(j) .AND. m == l) expected_d = pot%dion(i, j)
82 2736 : IF (ABS(projectors%d(p, q) - expected_d) > 1.e-15_dp) ERROR STOP "UPF D expansion failed"
83 : END DO
84 : END DO
85 : END DO
86 : END DO
87 6 : ld = MAX(na*ncoset(la), nb*ncoset(lb))
88 6 : ldai = ncoset(la + nder)
89 0 : ALLOCATE (sab(ld, ld*ncoset(nder)), work(ldai, ldai, ncoset(nder)), &
90 0 : moments(na*ncoset(la), nb*ncoset(lb), ncoset(nder) - 1), &
91 0 : integrals(na*ncoset(la), projectors%nproj, ncoset(nder)), &
92 126 : reference(na*ncoset(la), projectors%nproj, ncoset(nder)))
93 30 : DO geometry = 1, 4
94 6 : SELECT CASE (geometry)
95 : CASE (1)
96 6 : rac = 0.0_dp
97 : CASE (2)
98 6 : rac = [0.2_dp, -0.5_dp, 0.7_dp]
99 : CASE (3)
100 6 : rac = [3.0_dp, 2.0_dp, -1.0_dp]
101 : CASE (4)
102 24 : rac = [-7.0_dp, 1.0_dp, 2.0_dp]
103 : END SELECT
104 24 : sab = 0.0_dp
105 : CALL overlap(la, 0, na, SPREAD(100.0_dp, 1, na), alpha, &
106 : lb, 0, nb, SPREAD(100.0_dp, 1, nb), eta, rac, NORM2(rac), &
107 96 : sab, nder, .TRUE., work, ldai)
108 264 : DO id = 1, ncoset(nder)
109 240 : p = 0
110 1464 : DO i = 1, nb
111 1200 : l = pot%lbeta(i)
112 1200 : first = (id - 1)*ld + (i - 1)*ncoset(lb) + ncoset(l - 1) + 1
113 : reference(:, p + 1:p + 2*l + 1, id) = &
114 15830880 : MATMUL(sab(:na*ncoset(la), first:first + nco(l) - 1), TRANSPOSE(orbtramat(l)%slm))
115 1440 : p = p + 2*l + 1
116 : END DO
117 : END DO
118 24 : CALL upf_projector_overlap(projectors, la, alpha, rac, nder, integrals)
119 264 : DO id = 1, ncoset(nder)
120 802800 : scale = MAX(1.0_dp, MAXVAL(ABS(reference(:, :, id))))
121 802800 : error = MAXVAL(ABS(integrals(:, :, id) - reference(:, :, id)))/scale
122 240 : IF (io_unit > 0) WRITE (io_unit, '(A,3I4,ES16.7)') &
123 120 : "UPF radial Gaussian overlap/derivative:", grid, geometry, id, error
124 264 : IF (error > 1.e-10_dp) ERROR STOP "Numerical UPF overlaps disagree with analytic Gaussian integrals"
125 : END DO
126 : ! Independent moment implementation uses the OS recurrence with a fixed reference.
127 24 : origin = [0.4_dp, 0.6_dp, -0.3_dp]
128 : CALL moment(la, na, alpha, SPREAD(100.0_dp, 1, na), 0, &
129 96 : lb, nb, eta, SPREAD(100.0_dp, 1, nb), nder, origin, origin + rac, moments)
130 240 : DO id = 2, ncoset(nder)
131 216 : p = 0
132 1320 : DO i = 1, nb
133 1080 : l = pot%lbeta(i)
134 1080 : shell = (i - 1)*ncoset(lb) + ncoset(l - 1) + 1
135 : reference(:, p + 1:p + 2*l + 1, id) = &
136 14247792 : MATMUL(moments(:, shell:shell + nco(l) - 1, id - 1), TRANSPOSE(orbtramat(l)%slm))
137 1296 : p = p + 2*l + 1
138 : END DO
139 : END DO
140 24 : CALL upf_projector_overlap(projectors, la, alpha, rac, nder, integrals, moment_origin=origin)
141 802824 : scale = MAX(1.0_dp, MAXVAL(ABS(reference)))
142 802824 : error = MAXVAL(ABS(integrals - reference))/scale
143 24 : IF (io_unit > 0) WRITE (io_unit, '(A,2I4,ES16.7)') "UPF radial Gaussian moments:", grid, geometry, error
144 30 : IF (error > 1.e-10_dp) ERROR STOP "Numerical UPF moments disagree with analytic Gaussian integrals"
145 : END DO
146 8 : DEALLOCATE (pot%r, pot%rab, pot%beta, sab, work, moments, integrals, reference)
147 : END DO
148 2 : CALL atom_release_upf(pot)
149 : ! Optional library audit checks real radial meshes, angular channels, and
150 : ! successive cache replacement. Finiteness is not a physical accuracy test.
151 4 : DO i = 1, COMMAND_ARGUMENT_COUNT()
152 2 : CALL GET_COMMAND_ARGUMENT(i, filename)
153 : ! do_regtest.py passes the CP2K source root before any optional UPF audit files.
154 2 : IF (i == 1) THEN
155 2 : INQUIRE (FILE=TRIM(filename)//"/tests/UNIT_TESTS", EXIST=source_root_argument)
156 2 : IF (source_root_argument) CYCLE
157 : END IF
158 0 : CALL atom_read_upf(pot, TRIM(filename))
159 0 : IF (pot%number_of_proj > 0) THEN
160 0 : CALL upf_projector_init(pot, la + nder, projectors)
161 0 : IF (.NOT. ALL(ieee_is_finite(projectors%d))) ERROR STOP "Non-finite original UPF D matrix"
162 0 : ALLOCATE (integrals(na*ncoset(la), projectors%nproj, ncoset(nder)))
163 0 : DO geometry = 1, 2
164 0 : rac = 0.0_dp
165 0 : IF (geometry == 2) rac = [0.3_dp, -0.5_dp, 1.0_dp]
166 0 : CALL upf_projector_overlap(projectors, la, alpha, rac, nder, integrals)
167 0 : IF (.NOT. ALL(ieee_is_finite(integrals))) ERROR STOP "Non-finite native UPF projector integrals"
168 : END DO
169 0 : IF (io_unit > 0) WRITE (io_unit, '(A,1X,A,1X,I6,1X,ES16.7)') "UPF native integrals:", TRIM(filename), &
170 0 : projectors%nproj, MAXVAL(ABS(integrals))
171 0 : DEALLOCATE (integrals)
172 : END IF
173 4 : CALL atom_release_upf(pot)
174 : END DO
175 2 : CALL cp2k_finalize()
176 2 : END PROGRAM upf_projector_integrals_unittest
|
