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 : ! **************************************************************************************************
9 : !> \brief Common framework for using eigenvectors of a Fock matrix as PAO basis.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_param_fock
13 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
14 : dbcsr_get_info
15 : USE kinds, ONLY: dp
16 : USE mathlib, ONLY: diamat_all
17 : USE pao_types, ONLY: pao_env_type
18 : #include "./base/base_uses.f90"
19 :
20 : IMPLICIT NONE
21 :
22 : PRIVATE
23 :
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_fock'
25 :
26 : PUBLIC :: pao_calc_U_block_fock
27 :
28 : CONTAINS
29 :
30 : ! **************************************************************************************************
31 : !> \brief Calculate new matrix U and optinally its gradient G
32 : !> \param pao ...
33 : !> \param iatom ...
34 : !> \param V ...
35 : !> \param U ...
36 : !> \param penalty ...
37 : !> \param gap ...
38 : !> \param evals ...
39 : !> \param M1 ...
40 : !> \param G ...
41 : ! **************************************************************************************************
42 15505 : SUBROUTINE pao_calc_U_block_fock(pao, iatom, V, U, penalty, gap, evals, M1, G)
43 : TYPE(pao_env_type), POINTER :: pao
44 : INTEGER, INTENT(IN) :: iatom
45 : REAL(dp), DIMENSION(:, :), POINTER :: V, U
46 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
47 : REAL(dp), INTENT(OUT) :: gap
48 : REAL(dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: evals
49 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: M1, G
50 :
51 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_U_block_fock'
52 :
53 : INTEGER :: handle, i, j, m, n
54 15505 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pao, blk_sizes_pri
55 : LOGICAL :: found
56 : REAL(dp) :: alpha, beta, denom, diff
57 15505 : REAL(dp), DIMENSION(:), POINTER :: H_evals
58 15505 : REAL(dp), DIMENSION(:, :), POINTER :: block_N, D1, D2, H, H0, H_evecs, M2, M3, &
59 15505 : M4, M5
60 :
61 15505 : CALL timeset(routineN, handle)
62 :
63 15505 : CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=H0, found=found)
64 15505 : CPASSERT(ASSOCIATED(H0))
65 15505 : CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_N, found=found)
66 15505 : CPASSERT(ASSOCIATED(block_N))
67 1019353 : IF (MAXVAL(ABS(V - TRANSPOSE(V))) > 1e-14_dp) CPABORT("Expect symmetric matrix")
68 :
69 : ! figure out basis sizes
70 15505 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
71 15505 : n = blk_sizes_pri(iatom) ! size of primary basis
72 15505 : m = blk_sizes_pao(iatom) ! size of pao basis
73 :
74 : ! calculate H in the orthonormal basis
75 62020 : ALLOCATE (H(n, n))
76 60798414 : H = MATMUL(MATMUL(block_N, H0 + V), block_N)
77 :
78 : ! diagonalize H
79 77525 : ALLOCATE (H_evals(n), H_evecs(n, n))
80 2023201 : H_evecs = H
81 15505 : CALL diamat_all(H_evecs, H_evals)
82 :
83 : ! the eigenvectors of H become the rotation matrix U
84 2023201 : U = H_evecs
85 :
86 : ! copy eigenvectors around the gap from H_evals into evals array
87 15505 : IF (PRESENT(evals)) THEN
88 12771 : CPASSERT(MOD(SIZE(evals), 2) == 0) ! gap will be exactely in the middle
89 12771 : i = SIZE(evals)/2
90 12771 : j = MIN(m, i)
91 53340 : evals(1 + i - j:i) = H_evals(1 + m - j:m) ! eigenvalues below gap
92 12771 : j = MIN(n - m, i)
93 58643 : evals(i:i + j) = H_evals(m:m + j) ! eigenvalues above gap
94 : END IF
95 :
96 : ! calculate homo-lumo gap (it's useful for detecting numerical issues)
97 15505 : gap = HUGE(dp)
98 15505 : IF (m < n) THEN
99 : ! catch special case n==m
100 14958 : gap = H_evals(m + 1) - H_evals(m)
101 : END IF
102 :
103 15505 : IF (PRESENT(penalty)) THEN
104 : ! penalty terms: occupied and virtual eigenvalues repel each other
105 15185 : alpha = pao%penalty_strength
106 15185 : beta = pao%penalty_dist
107 67624 : DO i = 1, m
108 222208 : DO j = m + 1, n
109 154584 : diff = H_evals(i) - H_evals(j)
110 207023 : penalty = penalty + alpha*EXP(-(diff/beta)**2)
111 : END DO
112 : END DO
113 :
114 : ! regularization energy
115 1007619 : penalty = penalty + pao%regularization*SUM(V**2)
116 : END IF
117 :
118 15505 : IF (PRESENT(G)) THEN ! TURNING POINT (if calc grad) -------------------------
119 :
120 2465 : CPASSERT(PRESENT(M1))
121 :
122 : ! calculate derivatives between eigenvectors of H
123 22185 : ALLOCATE (D1(n, n), M2(n, n), M3(n, n), M4(n, n))
124 19024 : DO i = 1, n
125 156847 : DO j = 1, n
126 : ! ignore changes among occupied or virtual eigenvectors
127 : ! They will get filtered out by M2*D1 anyways, however this early
128 : ! intervention might stabilize numerics in the case of level-crossings.
129 154382 : IF (i <= m .EQV. j <= m) THEN
130 84603 : D1(i, j) = 0.0_dp
131 : ELSE
132 53220 : denom = H_evals(i) - H_evals(j)
133 53220 : IF (ABS(denom) > 1e-9_dp) THEN ! avoid division by zero
134 53220 : D1(i, j) = 1.0_dp/denom
135 : ELSE
136 0 : D1(i, j) = SIGN(1e+9_dp, denom)
137 : END IF
138 : END IF
139 : END DO
140 : END DO
141 2465 : IF (ASSOCIATED(M1)) THEN
142 3431248 : M2 = MATMUL(TRANSPOSE(M1), H_evecs)
143 : ELSE
144 0 : M2 = 0.0_dp
145 : END IF
146 311229 : M3 = M2*D1 ! Hadamard product
147 6244968 : M4 = MATMUL(MATMUL(H_evecs, M3), TRANSPOSE(H_evecs))
148 :
149 : ! gradient contribution from penalty terms
150 2465 : IF (PRESENT(penalty)) THEN
151 7293 : ALLOCATE (D2(n, n))
152 155769 : D2 = 0.0_dp
153 18818 : DO i = 1, n
154 155769 : DO j = 1, n
155 136951 : IF (i <= m .EQV. j <= m) CYCLE
156 52944 : diff = H_evals(i) - H_evals(j)
157 153338 : D2(i, i) = D2(i, i) - 2.0_dp*alpha*diff/beta**2*EXP(-(diff/beta)**2)
158 : END DO
159 : END DO
160 8876809 : M4 = M4 + MATMUL(MATMUL(H_evecs, D2), TRANSPOSE(H_evecs))
161 2431 : DEALLOCATE (D2)
162 : END IF
163 :
164 : ! dH / dV, return to non-orthonormal basis
165 7395 : ALLOCATE (M5(n, n))
166 11867478 : M5 = MATMUL(MATMUL(block_N, M4), block_N)
167 :
168 : ! add regularization gradient
169 2465 : IF (PRESENT(penalty)) THEN
170 309107 : M5 = M5 + 2.0_dp*pao%regularization*V
171 : END IF
172 :
173 : ! symmetrize
174 311229 : G = 0.5_dp*(M5 + TRANSPOSE(M5)) ! the final gradient
175 :
176 2465 : DEALLOCATE (D1, M2, M3, M4, M5)
177 : END IF
178 :
179 15505 : DEALLOCATE (H, H_evals, H_evecs)
180 :
181 15505 : CALL timestop(handle)
182 31010 : END SUBROUTINE pao_calc_U_block_fock
183 :
184 25331 : END MODULE pao_param_fock
|
