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 2 : PROGRAM kpsym_unittest
8 2 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
9 : USE cryssym, ONLY: crys_sym_gen,&
10 : csym_type,&
11 : kpoint_gen,&
12 : release_csym_type
13 : USE kinds, ONLY: dp
14 : USE libcp2k, ONLY: cp2k_finalize,&
15 : cp2k_init
16 : USE physcon, ONLY: bohr
17 :
18 : IMPLICIT NONE
19 : REAL(dp), PARAMETER :: ice_ix(3, 36) = RESHAPE([ &
20 : -0.006830_dp, 0.340009_dp, 0.190161_dp, &
21 : 0.113981_dp, 0.164412_dp, 0.283330_dp, &
22 : 0.309270_dp, 0.370346_dp, 0.106284_dp, &
23 : 0.159991_dp, 0.493170_dp, 0.440161_dp, &
24 : 0.335588_dp, 0.613981_dp, 0.533330_dp, &
25 : 0.129654_dp, 0.809270_dp, 0.356284_dp, &
26 : 0.840009_dp, 0.506830_dp, -0.059839_dp, &
27 : 0.664412_dp, 0.386019_dp, 0.033330_dp, &
28 : 0.870346_dp, 0.190730_dp, 0.856284_dp, &
29 : 0.506830_dp, 0.840009_dp, 1.059839_dp, &
30 : 0.386019_dp, 0.664412_dp, 0.966670_dp, &
31 : 0.190730_dp, 0.870346_dp, 0.143716_dp, &
32 : 0.493170_dp, 0.159991_dp, 0.559839_dp, &
33 : 0.613981_dp, 0.335588_dp, 0.466670_dp, &
34 : 0.809270_dp, 0.129654_dp, 0.643716_dp, &
35 : 0.370346_dp, 0.309270_dp, -0.106284_dp, &
36 : 0.659991_dp, 1.006830_dp, 0.309839_dp, &
37 : 0.835588_dp, 0.886019_dp, 0.216670_dp, &
38 : 0.629654_dp, 0.690730_dp, 0.393716_dp, &
39 : 1.006830_dp, 0.659991_dp, 0.690161_dp, &
40 : 0.886019_dp, 0.835588_dp, 0.783330_dp, &
41 : 0.690730_dp, 0.629654_dp, 0.606284_dp, &
42 : 0.340009_dp, -0.006830_dp, 0.809839_dp, &
43 : 0.164412_dp, 0.113981_dp, 0.716670_dp, &
44 : 0.114088_dp, 0.310030_dp, 0.267366_dp, &
45 : 0.189970_dp, 0.614088_dp, 0.517366_dp, &
46 : 0.810030_dp, 0.385912_dp, 0.017365_dp, &
47 : 0.385912_dp, 0.810030_dp, 0.982634_dp, &
48 : 0.614088_dp, 0.189970_dp, 0.482634_dp, &
49 : 0.689970_dp, 0.885912_dp, 0.232634_dp, &
50 : 0.885912_dp, 0.689970_dp, 0.767366_dp, &
51 : 0.310030_dp, 0.114088_dp, 0.732634_dp, &
52 : 0.404256_dp, 0.404256_dp, 0.000000_dp, &
53 : 0.095744_dp, 0.904256_dp, 0.250000_dp, &
54 : 0.904256_dp, 0.095744_dp, 0.750000_dp, &
55 : 0.595744_dp, 0.595744_dp, 0.500000_dp], [3, 36])
56 : REAL(dp) :: cell(3, 3), coords(3, 36), shifted(3, 36)
57 : INTEGER :: i, io_unit
58 : INTEGER :: types(36)
59 :
60 2 : CALL cp2k_init()
61 2 : io_unit = cp_logger_get_default_io_unit()
62 2 : cell = 0.0_dp
63 2 : cell(1, 1) = 6.784903_dp*bohr
64 2 : cell(2, 2) = 6.784903_dp*bohr
65 2 : cell(3, 3) = 6.807916_dp*bohr
66 50 : types(:24) = 1
67 26 : types(25:) = 2
68 2 : coords = ice_ix
69 :
70 : ! Rounded coordinates put seven operations within a few ulps of the boundary.
71 2 : CALL check_mesh(coords, cell, 1.0e-6_dp, 6, .FALSE., .TRUE.)
72 2 : CALL check_mesh(coords, cell, 1.0e-7_dp, 14, .FALSE., .TRUE.)
73 290 : shifted = coords + SPREAD([0.137_dp, -0.231_dp, 0.319_dp], 2, 36)
74 2 : CALL check_mesh(shifted, cell, 1.0e-6_dp, 6, .FALSE., .TRUE.)
75 26 : CALL check_mesh(coords, 1.3_dp*cell, 1.0e-6_dp, 6, .FALSE., .TRUE.)
76 50 : DO i = 1, 24
77 194 : shifted(:, i) = coords(:, 25 - i)
78 : END DO
79 26 : DO i = 25, 36
80 98 : shifted(:, i) = coords(:, 61 - i)
81 : END DO
82 2 : CALL check_mesh(shifted, cell, 1.0e-6_dp, 6, .FALSE., .TRUE.)
83 :
84 : ! Genuine tolerance-induced nonclosure must still trigger the orbit guard.
85 2 : coords = ice_ix
86 2 : coords(3, 27) = coords(3, 27) + 1.0e-6_dp
87 2 : coords(3, 33) = coords(3, 33) + 0.75e-6_dp
88 2 : coords(3, 36) = coords(3, 36) - 0.75e-6_dp
89 2 : CALL check_mesh(coords, cell, 1.0e-6_dp, 14, .TRUE., .FALSE.)
90 2 : IF (io_unit > 0) WRITE (io_unit, *) "K290 tolerance-boundary and orbit-guard tests passed"
91 2 : CALL cp2k_finalize()
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief Check mesh reduction and atomic operation consistency without an SCF calculation.
97 : !> \param x fractional atomic coordinates
98 : !> \param h cell matrix in bohr
99 : !> \param delta K290 symmetry tolerance
100 : !> \param expected expected irreducible k-point count
101 : !> \param fallback expected inversion fallback
102 : !> \param closed whether to check closure of the retained atomic operations
103 : ! **************************************************************************************************
104 12 : SUBROUTINE check_mesh(x, h, delta, expected, fallback, closed)
105 : REAL(dp), INTENT(IN) :: x(:, :), h(3, 3), delta
106 : INTEGER, INTENT(IN) :: expected
107 : LOGICAL, INTENT(IN) :: fallback, closed
108 :
109 : INTEGER :: a, i, j, k
110 : LOGICAL :: found
111 : REAL(dp) :: product(3, 3)
112 252 : TYPE(csym_type) :: csym
113 :
114 12 : CALL crys_sym_gen(csym, x, types, h, delta=delta, use_spglib=.FALSE.)
115 12 : CALL kpoint_gen(csym, [3, 3, 3], symm=.TRUE., shift=[0.0_dp, 0.0_dp, 0.0_dp])
116 12 : IF (csym%nkpoint /= expected) ERROR STOP "Incorrect irreducible mesh size"
117 12 : IF (csym%inversion_only .NEQV. fallback) ERROR STOP "Incorrect orbit-guard decision"
118 116 : IF (ABS(SUM(csym%wkpoint) - 27.0_dp) > 1.0e-12_dp) ERROR STOP "Incorrect weight sum"
119 984 : IF (ANY(csym%kplink <= 0)) ERROR STOP "Incomplete mesh mapping"
120 90 : DO k = 1, SIZE(csym%rt, 3)
121 2886 : DO a = 1, SIZE(types)
122 103974 : IF (COUNT(csym%f0(:, k) == a) /= 1) ERROR STOP "Nonbijective atom mapping"
123 : END DO
124 5784 : IF (ANY(types(csym%f0(:, k)) /= types)) ERROR STOP "Atom types not preserved"
125 : END DO
126 12 : IF (closed) THEN
127 76 : DO i = 1, csym%nrtot
128 590 : DO j = 1, csym%nrtot
129 34438 : product = MATMUL(csym%rt(:, :, i), csym%rt(:, :, j))
130 514 : found = .FALSE.
131 4612 : DO k = 1, csym%nrtot
132 53274 : IF (MAXVAL(ABS(product - csym%rt(:, :, k))) > 1.0e-12_dp) CYCLE
133 37522 : IF (ANY(csym%f0(csym%f0(:, j), i) /= csym%f0(:, k))) CYCLE
134 4612 : found = .TRUE.
135 : END DO
136 580 : IF (.NOT. found) ERROR STOP "Operation group not closed"
137 : END DO
138 : END DO
139 : END IF
140 12 : CALL release_csym_type(csym)
141 14 : END SUBROUTINE check_mesh
142 : END PROGRAM kpsym_unittest
|
