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 Original matrix exponential parametrization
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_param_exp
13 : USE basis_set_types, ONLY: gto_basis_set_type
14 : USE cp_dbcsr_api, ONLY: &
15 : dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
16 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
17 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
18 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_diag_blocks
19 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
20 : USE kinds, ONLY: dp
21 : USE mathlib, ONLY: diag_antisym,&
22 : diamat_all
23 : USE pao_param_methods, ONLY: pao_calc_AB_from_U,&
24 : pao_calc_grad_lnv_wrt_U
25 : USE pao_potentials, ONLY: pao_guess_initial_potential
26 : USE pao_types, ONLY: pao_env_type
27 : USE qs_environment_types, ONLY: get_qs_env,&
28 : qs_environment_type
29 : USE qs_kind_types, ONLY: get_qs_kind,&
30 : qs_kind_type
31 : #include "./base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_exp'
38 :
39 : PUBLIC :: pao_param_init_exp, pao_param_finalize_exp, pao_calc_AB_exp
40 : PUBLIC :: pao_param_count_exp, pao_param_initguess_exp
41 :
42 : CONTAINS
43 :
44 : ! **************************************************************************************************
45 : !> \brief Initialize matrix exponential parametrization
46 : !> \param pao ...
47 : !> \param qs_env ...
48 : ! **************************************************************************************************
49 24 : SUBROUTINE pao_param_init_exp(pao, qs_env)
50 : TYPE(pao_env_type), POINTER :: pao
51 : TYPE(qs_environment_type), POINTER :: qs_env
52 :
53 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_init_exp'
54 :
55 : INTEGER :: acol, arow, handle, iatom, N
56 : LOGICAL :: found
57 24 : REAL(dp), DIMENSION(:), POINTER :: H_evals
58 24 : REAL(dp), DIMENSION(:, :), POINTER :: block_H, block_H0, block_N, block_U0, &
59 24 : block_V0, H_evecs
60 : TYPE(dbcsr_iterator_type) :: iter
61 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
62 :
63 24 : CALL timeset(routineN, handle)
64 :
65 24 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
66 :
67 : ! allocate matrix_U0
68 : CALL dbcsr_create(pao%matrix_U0, &
69 : name="PAO matrix_U0", &
70 : matrix_type="N", &
71 : dist=pao%diag_distribution, &
72 24 : template=matrix_s(1)%matrix)
73 24 : CALL dbcsr_reserve_diag_blocks(pao%matrix_U0)
74 :
75 : ! diagonalize each block of H0 and store eigenvectors in U0
76 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env) &
77 24 : !$OMP PRIVATE(iter,arow,acol,iatom,N,found,block_H0,block_V0,block_N,block_H,block_U0,H_evecs,H_evals)
78 : CALL dbcsr_iterator_start(iter, pao%matrix_U0)
79 : DO WHILE (dbcsr_iterator_blocks_left(iter))
80 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_U0)
81 : iatom = arow; CPASSERT(arow == acol)
82 : CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=block_H0, found=found)
83 : CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_N, found=found)
84 : CPASSERT(ASSOCIATED(block_H0) .AND. ASSOCIATED(block_N))
85 : N = SIZE(block_U0, 1)
86 :
87 : ALLOCATE (block_V0(N, N))
88 : CALL pao_guess_initial_potential(qs_env, iatom, block_V0)
89 :
90 : ! construct H
91 : ALLOCATE (block_H(N, N))
92 : block_H = MATMUL(MATMUL(block_N, block_H0 + block_V0), block_N) ! transform into orthonormal basis
93 :
94 : ! diagonalize H
95 : ALLOCATE (H_evecs(N, N), H_evals(N))
96 : H_evecs = block_H
97 : CALL diamat_all(H_evecs, H_evals)
98 :
99 : ! use eigenvectors as initial guess
100 : block_U0 = H_evecs
101 :
102 : DEALLOCATE (block_H, H_evecs, H_evals, block_V0)
103 : END DO
104 : CALL dbcsr_iterator_stop(iter)
105 : !$OMP END PARALLEL
106 :
107 24 : IF (pao%precondition) THEN
108 0 : CPABORT("PAO preconditioning not supported for selected parametrization.")
109 : END IF
110 :
111 24 : CALL timestop(handle)
112 24 : END SUBROUTINE pao_param_init_exp
113 :
114 : ! **************************************************************************************************
115 : !> \brief Finalize exponential parametrization
116 : !> \param pao ...
117 : ! **************************************************************************************************
118 24 : SUBROUTINE pao_param_finalize_exp(pao)
119 : TYPE(pao_env_type), POINTER :: pao
120 :
121 24 : CALL dbcsr_release(pao%matrix_U0)
122 :
123 24 : END SUBROUTINE pao_param_finalize_exp
124 :
125 : ! **************************************************************************************************
126 : !> \brief Returns the number of parameters for given atomic kind
127 : !> \param qs_env ...
128 : !> \param ikind ...
129 : !> \param nparams ...
130 : ! **************************************************************************************************
131 128 : SUBROUTINE pao_param_count_exp(qs_env, ikind, nparams)
132 : TYPE(qs_environment_type), POINTER :: qs_env
133 : INTEGER, INTENT(IN) :: ikind
134 : INTEGER, INTENT(OUT) :: nparams
135 :
136 : INTEGER :: cols, pao_basis_size, pri_basis_size, &
137 : rows
138 : TYPE(gto_basis_set_type), POINTER :: basis_set
139 64 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
140 :
141 64 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
142 : CALL get_qs_kind(qs_kind_set(ikind), &
143 : basis_set=basis_set, &
144 64 : pao_basis_size=pao_basis_size)
145 64 : pri_basis_size = basis_set%nsgf
146 :
147 : ! we only consider rotations between occupied and virtuals
148 64 : rows = pao_basis_size
149 64 : cols = pri_basis_size - pao_basis_size
150 64 : nparams = rows*cols
151 :
152 64 : END SUBROUTINE pao_param_count_exp
153 :
154 : ! **************************************************************************************************
155 : !> \brief Fills matrix_X with an initial guess
156 : !> \param pao ...
157 : ! **************************************************************************************************
158 14 : SUBROUTINE pao_param_initguess_exp(pao)
159 : TYPE(pao_env_type), POINTER :: pao
160 :
161 14 : CALL dbcsr_set(pao%matrix_X, 0.0_dp) ! actual initial guess is matrix_U0
162 :
163 14 : END SUBROUTINE pao_param_initguess_exp
164 :
165 : ! **************************************************************************************************
166 : !> \brief Takes current matrix_X and calculates the matrices A and B.
167 : !> \param pao ...
168 : !> \param qs_env ...
169 : !> \param ls_scf_env ...
170 : !> \param gradient ...
171 : ! **************************************************************************************************
172 2710 : SUBROUTINE pao_calc_AB_exp(pao, qs_env, ls_scf_env, gradient)
173 : TYPE(pao_env_type), POINTER :: pao
174 : TYPE(qs_environment_type), POINTER :: qs_env
175 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
176 : LOGICAL, INTENT(IN) :: gradient
177 :
178 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_AB_exp'
179 :
180 : INTEGER :: handle
181 2710 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
182 : TYPE(dbcsr_type) :: matrix_M, matrix_U
183 :
184 2710 : CALL timeset(routineN, handle)
185 2710 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
186 2710 : CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
187 2710 : CALL dbcsr_reserve_diag_blocks(matrix_U)
188 :
189 : !TODO: move this condition into pao_calc_U, use matrix_N as template
190 2710 : IF (gradient) THEN
191 488 : CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
192 488 : CALL pao_calc_U_exp(pao, matrix_U, matrix_M, pao%matrix_G)
193 488 : CALL dbcsr_release(matrix_M)
194 : ELSE
195 2222 : CALL pao_calc_U_exp(pao, matrix_U)
196 : END IF
197 :
198 2710 : CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
199 2710 : CALL dbcsr_release(matrix_U)
200 2710 : CALL timestop(handle)
201 2710 : END SUBROUTINE pao_calc_AB_exp
202 :
203 : ! **************************************************************************************************
204 : !> \brief Calculate new matrix U and optionally its gradient G
205 : !> \param pao ...
206 : !> \param matrix_U ...
207 : !> \param matrix_M ...
208 : !> \param matrix_G ...
209 : ! **************************************************************************************************
210 2710 : SUBROUTINE pao_calc_U_exp(pao, matrix_U, matrix_M, matrix_G)
211 : TYPE(pao_env_type), POINTER :: pao
212 : TYPE(dbcsr_type) :: matrix_U
213 : TYPE(dbcsr_type), OPTIONAL :: matrix_M, matrix_G
214 :
215 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_U_exp'
216 :
217 : COMPLEX(dp) :: denom
218 2710 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
219 2710 : COMPLEX(dp), DIMENSION(:, :), POINTER :: block_D, evecs
220 : INTEGER :: acol, arow, handle, i, iatom, j, k, M, &
221 : N, nparams
222 2710 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pao, blk_sizes_pri
223 : LOGICAL :: found
224 2710 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_G_full, block_M, &
225 2710 : block_tmp, block_U, block_U0, block_X, &
226 2710 : block_X_full
227 : TYPE(dbcsr_iterator_type) :: iter
228 :
229 2710 : CALL timeset(routineN, handle)
230 :
231 2710 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
232 :
233 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_U,matrix_M,matrix_G,blk_sizes_pri,blk_sizes_pao) &
234 : !$OMP PRIVATE(iter,arow,acol,iatom,N,M,nparams,i,j,k,found) &
235 : !$OMP PRIVATE(block_X,block_U,block_U0,block_X_full,evals,evecs) &
236 2710 : !$OMP PRIVATE(block_M,block_G,block_D,block_tmp,block_G_full,denom)
237 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
238 : DO WHILE (dbcsr_iterator_blocks_left(iter))
239 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
240 : iatom = arow; CPASSERT(arow == acol)
241 : CALL dbcsr_get_block_p(matrix=matrix_U, row=iatom, col=iatom, block=block_U, found=found)
242 : CPASSERT(ASSOCIATED(block_U))
243 : CALL dbcsr_get_block_p(matrix=pao%matrix_U0, row=iatom, col=iatom, block=block_U0, found=found)
244 : CPASSERT(ASSOCIATED(block_U0))
245 :
246 : N = blk_sizes_pri(iatom) ! size of primary basis
247 : M = blk_sizes_pao(iatom) ! size of pao basis
248 : nparams = SIZE(block_X, 1)
249 :
250 : ! block_X stores only rotations between occupied and virtuals
251 : ! hence, we first have to build the full anti-symmetric exponent block
252 : ALLOCATE (block_X_full(N, N))
253 : block_X_full(:, :) = 0.0_dp
254 : DO i = 1, nparams
255 : block_X_full(MOD(i - 1, M) + 1, M + (i - 1)/M + 1) = +block_X(i, 1)
256 : block_X_full(M + (i - 1)/M + 1, MOD(i - 1, M) + 1) = -block_X(i, 1)
257 : END DO
258 :
259 : ! diagonalize block_X_full
260 : ALLOCATE (evals(N), evecs(N, N))
261 : CALL diag_antisym(block_X_full, evecs, evals)
262 :
263 : ! construct rotation matrix
264 : block_U(:, :) = 0.0_dp
265 : DO k = 1, N
266 : DO i = 1, N
267 : DO j = 1, N
268 : block_U(i, j) = block_U(i, j) + REAL(EXP(evals(k))*evecs(i, k)*CONJG(evecs(j, k)), dp)
269 : END DO
270 : END DO
271 : END DO
272 :
273 : block_U = MATMUL(block_U0, block_U) ! prepend initial guess rotation
274 :
275 : ! TURNING POINT (if calc grad) ------------------------------------------
276 : IF (PRESENT(matrix_G)) THEN
277 : CPASSERT(PRESENT(matrix_M))
278 :
279 : CALL dbcsr_get_block_p(matrix=pao%matrix_G, row=iatom, col=iatom, block=block_G, found=found)
280 : CPASSERT(ASSOCIATED(block_G))
281 : CALL dbcsr_get_block_p(matrix=matrix_M, row=iatom, col=iatom, block=block_M, found=found)
282 : ! don't check ASSOCIATED(block_M), it might have been filtered out.
283 :
284 : ALLOCATE (block_D(N, N), block_tmp(N, N), block_G_full(N, N))
285 : DO i = 1, N
286 : DO j = 1, N
287 : denom = evals(i) - evals(j)
288 : IF (i == j) THEN
289 : block_D(i, i) = EXP(evals(i)) ! diagonal elements
290 : ELSE IF (ABS(denom) > 1e-10_dp) THEN
291 : block_D(i, j) = (EXP(evals(i)) - EXP(evals(j)))/denom
292 : ELSE
293 : block_D(i, j) = 1.0_dp ! limit according to L'Hospital's rule
294 : END IF
295 : END DO
296 : END DO
297 :
298 : IF (ASSOCIATED(block_M)) THEN
299 : block_tmp = MATMUL(TRANSPOSE(block_U0), block_M)
300 : ELSE
301 : block_tmp = 0.0_dp
302 : END IF
303 : block_G_full = fold_derivatives(block_tmp, block_D, evecs)
304 :
305 : ! return only gradient for rotations between occupied and virtuals
306 : DO i = 1, nparams
307 : block_G(i, 1) = 2.0_dp*block_G_full(MOD(i - 1, M) + 1, M + (i - 1)/M + 1)
308 : END DO
309 :
310 : DEALLOCATE (block_D, block_tmp, block_G_full)
311 : END IF
312 :
313 : DEALLOCATE (block_X_full, evals, evecs)
314 :
315 : END DO
316 : CALL dbcsr_iterator_stop(iter)
317 : !$OMP END PARALLEL
318 :
319 2710 : CALL timestop(handle)
320 2710 : END SUBROUTINE pao_calc_U_exp
321 :
322 : ! **************************************************************************************************
323 : !> \brief Helper routine, for calculating derivatives
324 : !> \param M ...
325 : !> \param D ...
326 : !> \param R ...
327 : !> \return ...
328 : ! **************************************************************************************************
329 683 : FUNCTION fold_derivatives(M, D, R) RESULT(G)
330 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: M
331 : COMPLEX(dp), DIMENSION(:, :), INTENT(IN) :: D, R
332 : REAL(dp), DIMENSION(SIZE(M, 1), SIZE(M, 1)) :: G
333 :
334 683 : COMPLEX(dp), DIMENSION(:, :), POINTER :: F, RF, RM, RMR
335 : INTEGER :: n
336 683 : REAL(dp), DIMENSION(:, :), POINTER :: RFR
337 :
338 683 : n = SIZE(M, 1)
339 :
340 8879 : ALLOCATE (RM(n, n), RMR(n, n), F(n, n), RF(n, n), RFR(n, n))
341 :
342 641854 : RM = MATMUL(TRANSPOSE(CONJG(R)), TRANSPOSE(M))
343 1132635 : RMR = MATMUL(RM, R)
344 101626 : F = RMR*D !Hadamard product
345 641854 : RF = MATMUL(R, F)
346 1082505 : RFR = REAL(MATMUL(RF, TRANSPOSE(CONJG(R))), dp)
347 :
348 : ! gradient dE/dX has to be anti-symmetric
349 50813 : G = 0.5_dp*(TRANSPOSE(RFR) - RFR)
350 :
351 683 : DEALLOCATE (RM, RMR, F, RF, RFR)
352 683 : END FUNCTION fold_derivatives
353 :
354 683 : END MODULE pao_param_exp
|
