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 orbital_transformation_matrices_unittest
8 2 : USE kinds, ONLY: dp
9 : USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
10 : init_orbital_pointers
11 : USE orbital_transformation_matrices, ONLY: calculate_rotmat,&
12 : orbrotmat_type,&
13 : release_rotmat
14 : #include "../base/base_uses.f90"
15 :
16 : IMPLICIT NONE
17 :
18 : INTEGER, PARAMETER :: lmax = 5
19 : REAL(KIND=dp), PARAMETER :: eps = 1.0E-12_dp
20 : REAL(KIND=dp), DIMENSION(3, 3) :: rotmat
21 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrot => NULL()
22 :
23 2 : CALL init_orbital_pointers(lmax)
24 :
25 2 : rotmat = 0.0_dp
26 2 : rotmat(1, 1) = 1.0_dp
27 2 : rotmat(2, 2) = 1.0_dp
28 2 : rotmat(3, 3) = 1.0_dp
29 2 : CALL calculate_rotmat(orbrot, rotmat, lmax)
30 2 : CALL check_identity(orbrot)
31 2 : CALL check_orthogonal(orbrot)
32 :
33 2 : CALL make_z_rotation(0.7_dp, rotmat)
34 2 : CALL calculate_rotmat(orbrot, rotmat, lmax)
35 2 : CALL check_orthogonal(orbrot)
36 :
37 2 : CALL make_axis_rotation([1.0_dp, 2.0_dp, 3.0_dp], 0.37_dp, rotmat)
38 2 : CALL calculate_rotmat(orbrot, rotmat, lmax)
39 2 : CALL check_orthogonal(orbrot)
40 :
41 2 : CALL release_rotmat(orbrot)
42 2 : CALL deallocate_orbital_pointers()
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Check that the identity Cartesian rotation produces identity band matrices.
48 : !> \param orbrot ...
49 : ! **************************************************************************************************
50 2 : SUBROUTINE check_identity(orbrot)
51 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrot
52 :
53 : INTEGER :: i, l, n
54 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ref
55 :
56 18 : DO l = LBOUND(orbrot, 1), UBOUND(orbrot, 1)
57 12 : n = SIZE(orbrot(l)%mat, 1)
58 48 : ALLOCATE (ref(n, n))
59 12 : ref = 0.0_dp
60 84 : DO i = 1, n
61 84 : ref(i, i) = 1.0_dp
62 : END DO
63 656 : IF (MAXVAL(ABS(orbrot(l)%mat - ref)) > eps) THEN
64 0 : CPABORT("Bad identity rotation")
65 : END IF
66 14 : DEALLOCATE (ref)
67 : END DO
68 :
69 2 : END SUBROUTINE check_identity
70 :
71 : ! **************************************************************************************************
72 : !> \brief Check orthogonality of all band rotation matrices.
73 : !> \param orbrot ...
74 : ! **************************************************************************************************
75 6 : SUBROUTINE check_orthogonal(orbrot)
76 : TYPE(orbrotmat_type), DIMENSION(:), POINTER :: orbrot
77 :
78 : INTEGER :: i, j, k, l, n
79 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ident, prod
80 :
81 54 : DO l = LBOUND(orbrot, 1), UBOUND(orbrot, 1)
82 36 : n = SIZE(orbrot(l)%mat, 1)
83 216 : ALLOCATE (ident(n, n), prod(n, n))
84 36 : ident = 0.0_dp
85 252 : DO i = 1, n
86 252 : ident(i, i) = 1.0_dp
87 : END DO
88 36 : prod = 0.0_dp
89 252 : DO j = 1, n
90 1968 : DO k = 1, n
91 17268 : DO i = 1, n
92 17052 : prod(i, j) = prod(i, j) + orbrot(l)%mat(i, k)*orbrot(l)%mat(j, k)
93 : END DO
94 : END DO
95 : END DO
96 1968 : IF (MAXVAL(ABS(prod - ident)) > eps) THEN
97 0 : CPABORT("Bad rotation orthogonality")
98 : END IF
99 42 : DEALLOCATE (ident, prod)
100 : END DO
101 :
102 6 : END SUBROUTINE check_orthogonal
103 :
104 : ! **************************************************************************************************
105 : !> \brief Build a Cartesian rotation around the z axis.
106 : !> \param angle ...
107 : !> \param rotmat ...
108 : ! **************************************************************************************************
109 2 : SUBROUTINE make_z_rotation(angle, rotmat)
110 : REAL(KIND=dp), INTENT(IN) :: angle
111 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: rotmat
112 :
113 : REAL(KIND=dp) :: c, s
114 :
115 2 : c = COS(angle)
116 2 : s = SIN(angle)
117 2 : rotmat = 0.0_dp
118 2 : rotmat(1, 1) = c
119 2 : rotmat(1, 2) = -s
120 2 : rotmat(2, 1) = s
121 2 : rotmat(2, 2) = c
122 2 : rotmat(3, 3) = 1.0_dp
123 :
124 2 : END SUBROUTINE make_z_rotation
125 :
126 : ! **************************************************************************************************
127 : !> \brief Build a Cartesian rotation around an arbitrary axis.
128 : !> \param axis ...
129 : !> \param angle ...
130 : !> \param rotmat ...
131 : ! **************************************************************************************************
132 2 : SUBROUTINE make_axis_rotation(axis, angle, rotmat)
133 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: axis
134 : REAL(KIND=dp), INTENT(IN) :: angle
135 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: rotmat
136 :
137 : REAL(KIND=dp) :: c, norm, one_c, s, x, y, z
138 :
139 8 : norm = SQRT(SUM(axis**2))
140 2 : CPASSERT(norm > 0.0_dp)
141 2 : x = axis(1)/norm
142 2 : y = axis(2)/norm
143 2 : z = axis(3)/norm
144 2 : c = COS(angle)
145 2 : s = SIN(angle)
146 2 : one_c = 1.0_dp - c
147 :
148 2 : rotmat(1, 1) = c + x*x*one_c
149 2 : rotmat(1, 2) = x*y*one_c - z*s
150 2 : rotmat(1, 3) = x*z*one_c + y*s
151 2 : rotmat(2, 1) = y*x*one_c + z*s
152 2 : rotmat(2, 2) = c + y*y*one_c
153 2 : rotmat(2, 3) = y*z*one_c - x*s
154 2 : rotmat(3, 1) = z*x*one_c - y*s
155 2 : rotmat(3, 2) = z*y*one_c + x*s
156 2 : rotmat(3, 3) = c + z*z*one_c
157 :
158 2 : END SUBROUTINE make_axis_rotation
159 :
160 : END PROGRAM orbital_transformation_matrices_unittest
|
