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 4655 : 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 4655 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pao, blk_sizes_pri
55 : LOGICAL :: found
56 : REAL(dp) :: alpha, beta, denom, diff
57 4655 : REAL(dp), DIMENSION(:), POINTER :: H_evals
58 4655 : REAL(dp), DIMENSION(:, :), POINTER :: block_N, D1, D2, H, H0, H_evecs, M2, M3, &
59 4655 : M4, M5
60 :
61 4655 : CALL timeset(routineN, handle)
62 :
63 4655 : CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=H0, found=found)
64 4655 : CPASSERT(ASSOCIATED(H0))
65 4655 : CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_N, found=found)
66 4655 : CPASSERT(ASSOCIATED(block_N))
67 144969 : IF (MAXVAL(ABS(V - TRANSPOSE(V))) > 1e-14_dp) CPABORT("Expect symmetric matrix")
68 :
69 : ! figure out basis sizes
70 4655 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
71 4655 : n = blk_sizes_pri(iatom) ! size of primary basis
72 4655 : m = blk_sizes_pao(iatom) ! size of pao basis
73 :
74 : ! calculate H in the orthonormal basis
75 18620 : ALLOCATE (H(n, n))
76 4261994 : H = MATMUL(MATMUL(block_N, H0 + V), block_N)
77 :
78 : ! diagonalize H
79 23275 : ALLOCATE (H_evals(n), H_evecs(n, n))
80 285283 : H_evecs = H
81 4655 : CALL diamat_all(H_evecs, H_evals)
82 :
83 : ! the eigenvectors of H become the rotation matrix U
84 285283 : U = H_evecs
85 :
86 : ! copy eigenvectors around the gap from H_evals into evals array
87 4655 : IF (PRESENT(evals)) THEN
88 3670 : CPASSERT(MOD(SIZE(evals), 2) == 0) ! gap will be exactely in the middle
89 3670 : i = SIZE(evals)/2
90 3670 : j = MIN(m, i)
91 7348 : evals(1 + i - j:i) = H_evals(1 + m - j:m) ! eigenvalues below gap
92 3670 : j = MIN(n - m, i)
93 22042 : 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 4655 : gap = HUGE(dp)
98 4655 : IF (m < n) THEN
99 : ! catch special case n==m
100 4653 : gap = H_evals(m + 1) - H_evals(m)
101 : END IF
102 :
103 4655 : IF (PRESENT(penalty)) THEN
104 : ! penalty terms: occupied and virtual eigenvalues repel each other
105 4362 : alpha = pao%penalty_strength
106 4362 : beta = pao%penalty_dist
107 8724 : DO i = 1, m
108 26172 : DO j = m + 1, n
109 17448 : diff = H_evals(i) - H_evals(j)
110 21810 : penalty = penalty + alpha*EXP(-(diff/beta)**2)
111 : END DO
112 : END DO
113 :
114 : ! regularization energy
115 135222 : penalty = penalty + pao%regularization*SUM(V**2)
116 : END IF
117 :
118 4655 : IF (PRESENT(G)) THEN ! TURNING POINT (if calc grad) -------------------------
119 :
120 660 : CPASSERT(PRESENT(M1))
121 :
122 : ! calculate derivatives between eigenvectors of H
123 5940 : ALLOCATE (D1(n, n), M2(n, n), M3(n, n), M4(n, n))
124 3962 : DO i = 1, n
125 20484 : 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 19824 : IF (i <= m .EQV. j <= m) THEN
130 11238 : D1(i, j) = 0.0_dp
131 : ELSE
132 5284 : denom = H_evals(i) - H_evals(j)
133 5284 : IF (ABS(denom) > 1e-9_dp) THEN ! avoid division by zero
134 5284 : 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 660 : IF (ASSOCIATED(M1)) THEN
142 245980 : M2 = MATMUL(TRANSPOSE(M1), H_evecs)
143 : ELSE
144 0 : M2 = 0.0_dp
145 : END IF
146 40308 : M3 = M2*D1 ! Hadamard product
147 412664 : M4 = MATMUL(MATMUL(H_evecs, M3), TRANSPOSE(H_evecs))
148 :
149 : ! gradient contribution from penalty terms
150 660 : IF (PRESENT(penalty)) THEN
151 1878 : ALLOCATE (D2(n, n))
152 19406 : D2 = 0.0_dp
153 3756 : DO i = 1, n
154 19406 : DO j = 1, n
155 15650 : IF (i <= m .EQV. j <= m) CYCLE
156 5008 : diff = H_evals(i) - H_evals(j)
157 18780 : D2(i, i) = D2(i, i) - 2.0_dp*alpha*diff/beta**2*EXP(-(diff/beta)**2)
158 : END DO
159 : END DO
160 530222 : M4 = M4 + MATMUL(MATMUL(H_evecs, D2), TRANSPOSE(H_evecs))
161 626 : DEALLOCATE (D2)
162 : END IF
163 :
164 : ! dH / dV, return to non-orthonormal basis
165 1980 : ALLOCATE (M5(n, n))
166 744712 : M5 = MATMUL(MATMUL(block_N, M4), block_N)
167 :
168 : ! add regularization gradient
169 660 : IF (PRESENT(penalty)) THEN
170 38186 : M5 = M5 + 2.0_dp*pao%regularization*V
171 : END IF
172 :
173 : ! symmetrize
174 40308 : G = 0.5_dp*(M5 + TRANSPOSE(M5)) ! the final gradient
175 :
176 660 : DEALLOCATE (D1, M2, M3, M4, M5)
177 : END IF
178 :
179 4655 : DEALLOCATE (H, H_evals, H_evecs)
180 :
181 4655 : CALL timestop(handle)
182 9310 : END SUBROUTINE pao_calc_U_block_fock
183 :
184 7261 : END MODULE pao_param_fock
|
