Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 15493 : 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 15493 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pao, blk_sizes_pri
55 : LOGICAL :: found
56 : REAL(dp) :: alpha, beta, denom, diff
57 15493 : REAL(dp), DIMENSION(:), POINTER :: H_evals
58 15493 : REAL(dp), DIMENSION(:, :), POINTER :: block_N, D1, D2, H, H0, H_evecs, M2, M3, &
59 15493 : M4, M5
60 :
61 15493 : CALL timeset(routineN, handle)
62 :
63 15493 : CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=H0, found=found)
64 15493 : CPASSERT(ASSOCIATED(H0))
65 15493 : CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_N, found=found)
66 15493 : CPASSERT(ASSOCIATED(block_N))
67 1018373 : IF (MAXVAL(ABS(V - TRANSPOSE(V))) > 1e-14_dp) CPABORT("Expect symmetric matrix")
68 :
69 : ! figure out basis sizes
70 15493 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
71 15493 : n = blk_sizes_pri(iatom) ! size of primary basis
72 15493 : m = blk_sizes_pao(iatom) ! size of pao basis
73 :
74 : ! calculate H in the orthonormal basis
75 61972 : ALLOCATE (H(n, n))
76 61737654 : H = MATMUL(MATMUL(block_N, H0 + V), block_N)
77 :
78 : ! diagonalize H
79 77465 : ALLOCATE (H_evals(n), H_evecs(n, n))
80 2021253 : H_evecs = H
81 15493 : CALL diamat_all(H_evecs, H_evals)
82 :
83 : ! the eigenvectors of H become the rotation matrix U
84 2021253 : U = H_evecs
85 :
86 : ! copy eigenvectors around the gap from H_evals into evals array
87 15493 : IF (PRESENT(evals)) THEN
88 12759 : CPASSERT(MOD(SIZE(evals), 2) == 0) ! gap will be exactely in the middle
89 12759 : i = SIZE(evals)/2
90 12759 : j = MIN(m, i)
91 53276 : evals(1 + i - j:i) = H_evals(1 + m - j:m) ! eigenvalues below gap
92 12759 : j = MIN(n - m, i)
93 58611 : 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 15493 : gap = HUGE(dp)
98 15493 : IF (m < n) & ! catch special case n==m
99 14950 : gap = H_evals(m + 1) - H_evals(m)
100 :
101 15493 : IF (PRESENT(penalty)) THEN
102 : ! penalty terms: occupied and virtual eigenvalues repel each other
103 15173 : alpha = pao%penalty_strength
104 15173 : beta = pao%penalty_dist
105 67528 : DO i = 1, m
106 222080 : DO j = m + 1, n
107 154552 : diff = H_evals(i) - H_evals(j)
108 206907 : penalty = penalty + alpha*EXP(-(diff/beta)**2)
109 : END DO
110 : END DO
111 :
112 : ! regularization energy
113 1006639 : penalty = penalty + pao%regularization*SUM(V**2)
114 : END IF
115 :
116 15493 : IF (PRESENT(G)) THEN ! TURNING POINT (if calc grad) -------------------------
117 :
118 2465 : CPASSERT(PRESENT(M1))
119 :
120 : ! calculate derivatives between eigenvectors of H
121 22185 : ALLOCATE (D1(n, n), M2(n, n), M3(n, n), M4(n, n))
122 19024 : DO i = 1, n
123 156847 : DO j = 1, n
124 : ! ignore changes among occupied or virtual eigenvectors
125 : ! They will get filtered out by M2*D1 anyways, however this early
126 : ! intervention might stabilize numerics in the case of level-crossings.
127 154382 : IF (i <= m .EQV. j <= m) THEN
128 84603 : D1(i, j) = 0.0_dp
129 : ELSE
130 53220 : denom = H_evals(i) - H_evals(j)
131 53220 : IF (ABS(denom) > 1e-9_dp) THEN ! avoid division by zero
132 53220 : D1(i, j) = 1.0_dp/denom
133 : ELSE
134 0 : D1(i, j) = SIGN(1e+9_dp, denom)
135 : END IF
136 : END IF
137 : END DO
138 : END DO
139 2465 : IF (ASSOCIATED(M1)) THEN
140 3431248 : M2 = MATMUL(TRANSPOSE(M1), H_evecs)
141 : ELSE
142 0 : M2 = 0.0_dp
143 : END IF
144 311229 : M3 = M2*D1 ! Hadamard product
145 6399350 : M4 = MATMUL(MATMUL(H_evecs, M3), TRANSPOSE(H_evecs))
146 :
147 : ! gradient contribution from penalty terms
148 2465 : IF (PRESENT(penalty)) THEN
149 7293 : ALLOCATE (D2(n, n))
150 155769 : D2 = 0.0_dp
151 18818 : DO i = 1, n
152 155769 : DO j = 1, n
153 136951 : IF (i <= m .EQV. j <= m) CYCLE
154 52944 : diff = H_evals(i) - H_evals(j)
155 153338 : D2(i, i) = D2(i, i) - 2.0_dp*alpha*diff/beta**2*EXP(-(diff/beta)**2)
156 : END DO
157 : END DO
158 9183485 : M4 = M4 + MATMUL(MATMUL(H_evecs, D2), TRANSPOSE(H_evecs))
159 2431 : DEALLOCATE (D2)
160 : END IF
161 :
162 : ! dH / dV, return to non-orthonormal basis
163 7395 : ALLOCATE (M5(n, n))
164 12021860 : M5 = MATMUL(MATMUL(block_N, M4), block_N)
165 :
166 : ! add regularization gradient
167 2465 : IF (PRESENT(penalty)) &
168 309107 : M5 = M5 + 2.0_dp*pao%regularization*V
169 :
170 : ! symmetrize
171 311229 : G = 0.5_dp*(M5 + TRANSPOSE(M5)) ! the final gradient
172 :
173 2465 : DEALLOCATE (D1, M2, M3, M4, M5)
174 : END IF
175 :
176 15493 : DEALLOCATE (H, H_evals, H_evecs)
177 :
178 15493 : CALL timestop(handle)
179 30986 : END SUBROUTINE pao_calc_U_block_fock
180 :
181 25319 : END MODULE pao_param_fock
|
