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