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 a linear parametrization of the potential.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_param_linpot
13 : USE atomic_kind_types, ONLY: get_atomic_kind
14 : USE basis_set_types, ONLY: gto_basis_set_type
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
18 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
19 : dbcsr_p_type, dbcsr_release, dbcsr_type
20 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_diag_blocks
21 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
22 : USE kinds, ONLY: dp
23 : USE machine, ONLY: m_flush
24 : USE mathlib, ONLY: diamat_all
25 : USE message_passing, ONLY: mp_comm_type,&
26 : mp_para_env_type
27 : USE pao_input, ONLY: pao_fock_param,&
28 : pao_rotinv_param
29 : USE pao_linpot_full, ONLY: linpot_full_calc_terms,&
30 : linpot_full_count_terms
31 : USE pao_linpot_rotinv, ONLY: linpot_rotinv_calc_forces,&
32 : linpot_rotinv_calc_terms,&
33 : linpot_rotinv_count_terms
34 : USE pao_param_fock, ONLY: pao_calc_U_block_fock
35 : USE pao_param_methods, ONLY: pao_calc_AB_from_U,&
36 : pao_calc_grad_lnv_wrt_U
37 : USE pao_potentials, ONLY: pao_guess_initial_potential
38 : USE pao_types, ONLY: pao_env_type
39 : USE particle_types, ONLY: particle_type
40 : USE qs_environment_types, ONLY: get_qs_env,&
41 : qs_environment_type
42 : USE qs_kind_types, ONLY: get_qs_kind,&
43 : qs_kind_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 :
48 : PRIVATE
49 :
50 : PUBLIC :: pao_param_init_linpot, pao_param_finalize_linpot, pao_calc_AB_linpot
51 : PUBLIC :: pao_param_count_linpot, pao_param_initguess_linpot
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Initialize the linear potential parametrization
57 : !> \param pao ...
58 : !> \param qs_env ...
59 : ! **************************************************************************************************
60 234 : SUBROUTINE pao_param_init_linpot(pao, qs_env)
61 : TYPE(pao_env_type), POINTER :: pao
62 : TYPE(qs_environment_type), POINTER :: qs_env
63 :
64 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_init_linpot'
65 :
66 : INTEGER :: acol, arow, handle, iatom, ikind, N, &
67 : natoms, nterms
68 234 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, row_blk_size
69 234 : REAL(dp), DIMENSION(:, :), POINTER :: block_V_terms
70 234 : REAL(dp), DIMENSION(:, :, :), POINTER :: V_blocks
71 : TYPE(dbcsr_iterator_type) :: iter
72 : TYPE(dft_control_type), POINTER :: dft_control
73 : TYPE(mp_para_env_type), POINTER :: para_env
74 234 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
75 :
76 234 : CALL timeset(routineN, handle)
77 :
78 : CALL get_qs_env(qs_env, &
79 : para_env=para_env, &
80 : dft_control=dft_control, &
81 : particle_set=particle_set, &
82 234 : natom=natoms)
83 :
84 234 : IF (dft_control%nspins /= 1) CPABORT("open shell not yet implemented")
85 :
86 : ! figure out number of potential terms
87 936 : ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
88 714 : DO iatom = 1, natoms
89 480 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
90 480 : CALL pao_param_count_linpot(pao, qs_env, ikind, nterms)
91 714 : col_blk_size(iatom) = nterms
92 : END DO
93 :
94 : ! allocate matrix_V_terms
95 234 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri)
96 1428 : row_blk_size = blk_sizes_pri**2
97 : CALL dbcsr_create(pao%matrix_V_terms, &
98 : name="PAO matrix_V_terms", &
99 : dist=pao%diag_distribution, &
100 : matrix_type="N", &
101 : row_blk_size=row_blk_size, &
102 234 : col_blk_size=col_blk_size)
103 234 : CALL dbcsr_reserve_diag_blocks(pao%matrix_V_terms)
104 234 : DEALLOCATE (row_blk_size, col_blk_size)
105 :
106 : ! calculate, normalize, and store potential terms as rows of block_V_terms
107 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,blk_sizes_pri) &
108 234 : !$OMP PRIVATE(iter,arow,acol,iatom,N,nterms,block_V_terms,V_blocks)
109 : CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
110 : DO WHILE (dbcsr_iterator_blocks_left(iter))
111 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_V_terms)
112 : iatom = arow; CPASSERT(arow == acol)
113 : nterms = SIZE(block_V_terms, 2)
114 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
115 : N = blk_sizes_pri(iatom)
116 : CPASSERT(N*N == SIZE(block_V_terms, 1))
117 : ALLOCATE (V_blocks(N, N, nterms))
118 : CALL linpot_calc_terms(pao, qs_env, iatom, V_blocks)
119 : block_V_terms = RESHAPE(V_blocks, (/N*N, nterms/)) ! convert matrices into vectors
120 : DEALLOCATE (V_blocks)
121 : END DO
122 : CALL dbcsr_iterator_stop(iter)
123 : !$OMP END PARALLEL
124 :
125 234 : CALL pao_param_linpot_regularizer(pao)
126 :
127 234 : IF (pao%precondition) &
128 12 : CALL pao_param_linpot_preconditioner(pao)
129 :
130 234 : CALL para_env%sync() ! ensure that timestop is not called too early
131 :
132 234 : CALL timestop(handle)
133 234 : END SUBROUTINE pao_param_init_linpot
134 :
135 : ! **************************************************************************************************
136 : !> \brief Builds the regularization metric matrix_R
137 : !> \param pao ...
138 : ! **************************************************************************************************
139 234 : SUBROUTINE pao_param_linpot_regularizer(pao)
140 : TYPE(pao_env_type), POINTER :: pao
141 :
142 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_linpot_regularizer'
143 :
144 : INTEGER :: acol, arow, handle, i, iatom, j, k, &
145 : nterms
146 234 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
147 : LOGICAL :: found
148 : REAL(dp) :: v, w
149 234 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
150 234 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs
151 234 : REAL(dp), DIMENSION(:, :), POINTER :: block_R, V_terms
152 : TYPE(dbcsr_iterator_type) :: iter
153 :
154 234 : CALL timeset(routineN, handle)
155 :
156 234 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot regularizer"
157 :
158 234 : CALL dbcsr_get_info(pao%matrix_V_terms, col_blk_size=blk_sizes_nterms)
159 :
160 : ! build regularization metric
161 : CALL dbcsr_create(pao%matrix_R, &
162 : template=pao%matrix_V_terms, &
163 : matrix_type="N", &
164 : row_blk_size=blk_sizes_nterms, &
165 : col_blk_size=blk_sizes_nterms, &
166 234 : name="PAO matrix_R")
167 234 : CALL dbcsr_reserve_diag_blocks(pao%matrix_R)
168 :
169 : ! fill matrix_R
170 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
171 234 : !$OMP PRIVATE(iter,arow,acol,iatom,block_R,V_terms,found,nterms,S,S_evecs,S_evals,k,i,j,v,w)
172 : CALL dbcsr_iterator_start(iter, pao%matrix_R)
173 : DO WHILE (dbcsr_iterator_blocks_left(iter))
174 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_R)
175 : iatom = arow; CPASSERT(arow == acol)
176 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=V_terms, found=found)
177 : CPASSERT(ASSOCIATED(V_terms))
178 : nterms = SIZE(V_terms, 2)
179 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
180 :
181 : ! build overlap matrix
182 : ALLOCATE (S(nterms, nterms))
183 : S(:, :) = MATMUL(TRANSPOSE(V_terms), V_terms)
184 :
185 : ! diagonalize S
186 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
187 : S_evecs(:, :) = S
188 : CALL diamat_all(S_evecs, S_evals)
189 :
190 : block_R = 0.0_dp
191 : DO k = 1, nterms
192 : v = pao%linpot_regu_delta/S_evals(k)
193 : w = pao%linpot_regu_strength*MIN(1.0_dp, ABS(v))
194 : DO i = 1, nterms
195 : DO j = 1, nterms
196 : block_R(i, j) = block_R(i, j) + w*S_evecs(i, k)*S_evecs(j, k)
197 : END DO
198 : END DO
199 : END DO
200 :
201 : ! clean up
202 : DEALLOCATE (S, S_evals, S_evecs)
203 : END DO
204 : CALL dbcsr_iterator_stop(iter)
205 : !$OMP END PARALLEL
206 :
207 234 : CALL timestop(handle)
208 468 : END SUBROUTINE pao_param_linpot_regularizer
209 :
210 : ! **************************************************************************************************
211 : !> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
212 : !> \param pao ...
213 : ! **************************************************************************************************
214 12 : SUBROUTINE pao_param_linpot_preconditioner(pao)
215 : TYPE(pao_env_type), POINTER :: pao
216 :
217 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_linpot_preconditioner'
218 :
219 : INTEGER :: acol, arow, handle, i, iatom, j, k, &
220 : nterms
221 12 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
222 : LOGICAL :: found
223 : REAL(dp) :: eval_capped
224 12 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
225 12 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs
226 12 : REAL(dp), DIMENSION(:, :), POINTER :: block_precon, block_precon_inv, &
227 12 : block_V_terms
228 : TYPE(dbcsr_iterator_type) :: iter
229 :
230 12 : CALL timeset(routineN, handle)
231 :
232 12 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot preconditioner"
233 :
234 12 : CALL dbcsr_get_info(pao%matrix_V_terms, col_blk_size=blk_sizes_nterms)
235 :
236 : CALL dbcsr_create(pao%matrix_precon, &
237 : template=pao%matrix_V_terms, &
238 : matrix_type="N", &
239 : row_blk_size=blk_sizes_nterms, &
240 : col_blk_size=blk_sizes_nterms, &
241 12 : name="PAO matrix_precon")
242 12 : CALL dbcsr_reserve_diag_blocks(pao%matrix_precon)
243 :
244 12 : CALL dbcsr_create(pao%matrix_precon_inv, template=pao%matrix_precon, name="PAO matrix_precon_inv")
245 12 : CALL dbcsr_reserve_diag_blocks(pao%matrix_precon_inv)
246 :
247 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
248 12 : !$OMP PRIVATE(iter,arow,acol,iatom,block_V_terms,block_precon,block_precon_inv,found,nterms,S,S_evals,S_evecs,i,j,k,eval_capped)
249 : CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
250 : DO WHILE (dbcsr_iterator_blocks_left(iter))
251 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_V_terms)
252 : iatom = arow; CPASSERT(arow == acol)
253 : nterms = SIZE(block_V_terms, 2)
254 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
255 :
256 : CALL dbcsr_get_block_p(matrix=pao%matrix_precon, row=iatom, col=iatom, block=block_precon, found=found)
257 : CALL dbcsr_get_block_p(matrix=pao%matrix_precon_inv, row=iatom, col=iatom, block=block_precon_inv, found=found)
258 : CPASSERT(ASSOCIATED(block_precon))
259 : CPASSERT(ASSOCIATED(block_precon_inv))
260 :
261 : ALLOCATE (S(nterms, nterms))
262 : S(:, :) = MATMUL(TRANSPOSE(block_V_terms), block_V_terms)
263 :
264 : ! diagonalize S
265 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
266 : S_evecs(:, :) = S
267 : CALL diamat_all(S_evecs, S_evals)
268 :
269 : ! construct 1/Sqrt(S) and Sqrt(S)
270 : block_precon = 0.0_dp
271 : block_precon_inv = 0.0_dp
272 : DO k = 1, nterms
273 : eval_capped = MAX(pao%linpot_precon_delta, S_evals(k)) ! too small eigenvalues are hurtful
274 : DO i = 1, nterms
275 : DO j = 1, nterms
276 : block_precon(i, j) = block_precon(i, j) + S_evecs(i, k)*S_evecs(j, k)/SQRT(eval_capped)
277 : block_precon_inv(i, j) = block_precon_inv(i, j) + S_evecs(i, k)*S_evecs(j, k)*SQRT(eval_capped)
278 : END DO
279 : END DO
280 : END DO
281 :
282 : DEALLOCATE (S, S_evecs, S_evals)
283 : END DO
284 : CALL dbcsr_iterator_stop(iter)
285 : !$OMP END PARALLEL
286 :
287 12 : CALL timestop(handle)
288 24 : END SUBROUTINE pao_param_linpot_preconditioner
289 :
290 : ! **************************************************************************************************
291 : !> \brief Finalize the linear potential parametrization
292 : !> \param pao ...
293 : ! **************************************************************************************************
294 234 : SUBROUTINE pao_param_finalize_linpot(pao)
295 : TYPE(pao_env_type), POINTER :: pao
296 :
297 234 : CALL dbcsr_release(pao%matrix_V_terms)
298 234 : CALL dbcsr_release(pao%matrix_R)
299 :
300 234 : IF (pao%precondition) THEN
301 12 : CALL dbcsr_release(pao%matrix_precon)
302 12 : CALL dbcsr_release(pao%matrix_precon_inv)
303 : END IF
304 :
305 234 : END SUBROUTINE pao_param_finalize_linpot
306 :
307 : ! **************************************************************************************************
308 : !> \brief Returns the number of potential terms for given atomic kind
309 : !> \param pao ...
310 : !> \param qs_env ...
311 : !> \param ikind ...
312 : !> \param nparams ...
313 : ! **************************************************************************************************
314 1344 : SUBROUTINE pao_param_count_linpot(pao, qs_env, ikind, nparams)
315 : TYPE(pao_env_type), POINTER :: pao
316 : TYPE(qs_environment_type), POINTER :: qs_env
317 : INTEGER, INTENT(IN) :: ikind
318 : INTEGER, INTENT(OUT) :: nparams
319 :
320 : INTEGER :: pao_basis_size
321 : TYPE(gto_basis_set_type), POINTER :: basis_set
322 672 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
323 :
324 672 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
325 :
326 : CALL get_qs_kind(qs_kind_set(ikind), &
327 : basis_set=basis_set, &
328 672 : pao_basis_size=pao_basis_size)
329 :
330 672 : IF (pao_basis_size == basis_set%nsgf) THEN
331 26 : nparams = 0 ! pao disabled for iatom
332 :
333 : ELSE
334 754 : SELECT CASE (pao%parameterization)
335 : CASE (pao_fock_param)
336 646 : CALL linpot_full_count_terms(qs_env, ikind, nterms=nparams)
337 : CASE (pao_rotinv_param)
338 538 : CALL linpot_rotinv_count_terms(qs_env, ikind, nterms=nparams)
339 : CASE DEFAULT
340 646 : CPABORT("unknown parameterization")
341 : END SELECT
342 : END IF
343 :
344 672 : END SUBROUTINE pao_param_count_linpot
345 :
346 : ! **************************************************************************************************
347 : !> \brief Takes current matrix_X and calculates the matrices A and B.
348 : !> \param pao ...
349 : !> \param qs_env ...
350 : !> \param ls_scf_env ...
351 : !> \param gradient ...
352 : !> \param penalty ...
353 : !> \param forces ...
354 : ! **************************************************************************************************
355 8188 : SUBROUTINE pao_calc_AB_linpot(pao, qs_env, ls_scf_env, gradient, penalty, forces)
356 : TYPE(pao_env_type), POINTER :: pao
357 : TYPE(qs_environment_type), POINTER :: qs_env
358 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
359 : LOGICAL, INTENT(IN) :: gradient
360 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
361 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
362 :
363 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_AB_linpot'
364 :
365 : INTEGER :: handle
366 8188 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
367 : TYPE(dbcsr_type) :: matrix_M, matrix_U
368 :
369 8188 : CALL timeset(routineN, handle)
370 8188 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
371 8188 : CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
372 8188 : CALL dbcsr_reserve_diag_blocks(matrix_U)
373 :
374 : !TODO: move this condition into pao_calc_U, use matrix_N as template
375 8188 : IF (gradient) THEN
376 1616 : CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
377 3198 : CALL pao_calc_U_linpot(pao, qs_env, matrix_U, matrix_M, pao%matrix_G, penalty, forces)
378 1616 : CALL dbcsr_release(matrix_M)
379 : ELSE
380 6572 : CALL pao_calc_U_linpot(pao, qs_env, matrix_U, penalty=penalty)
381 : END IF
382 :
383 8188 : CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
384 8188 : CALL dbcsr_release(matrix_U)
385 8188 : CALL timestop(handle)
386 8188 : END SUBROUTINE pao_calc_AB_linpot
387 :
388 : ! **************************************************************************************************
389 : !> \brief Calculate new matrix U and optinally its gradient G
390 : !> \param pao ...
391 : !> \param qs_env ...
392 : !> \param matrix_U ...
393 : !> \param matrix_M ...
394 : !> \param matrix_G ...
395 : !> \param penalty ...
396 : !> \param forces ...
397 : ! **************************************************************************************************
398 8188 : SUBROUTINE pao_calc_U_linpot(pao, qs_env, matrix_U, matrix_M, matrix_G, penalty, forces)
399 : TYPE(pao_env_type), POINTER :: pao
400 : TYPE(qs_environment_type), POINTER :: qs_env
401 : TYPE(dbcsr_type) :: matrix_U
402 : TYPE(dbcsr_type), OPTIONAL :: matrix_M, matrix_G
403 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
404 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
405 :
406 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_U_linpot'
407 :
408 : INTEGER :: acol, arow, handle, iatom, kterm, n, &
409 : natoms, nterms
410 : LOGICAL :: found
411 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
412 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: evals
413 8188 : REAL(dp), DIMENSION(:), POINTER :: vec_M2, vec_V
414 8188 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_M1, block_M2, block_R, &
415 8188 : block_U, block_V, block_V_terms, &
416 8188 : block_X
417 8188 : REAL(dp), DIMENSION(:, :, :), POINTER :: M_blocks
418 : REAL(KIND=dp) :: regu_energy
419 : TYPE(dbcsr_iterator_type) :: iter
420 : TYPE(mp_comm_type) :: group
421 :
422 8188 : CALL timeset(routineN, handle)
423 :
424 8188 : CPASSERT(PRESENT(matrix_G) .EQV. PRESENT(matrix_M))
425 :
426 8188 : CALL get_qs_env(qs_env, natom=natoms)
427 40940 : ALLOCATE (gaps(natoms), evals(10, natoms)) ! printing 10 eigenvalues seems reasonable
428 244468 : evals(:, :) = 0.0_dp
429 29668 : gaps(:) = HUGE(1.0_dp)
430 8188 : regu_energy = 0.0_dp
431 8188 : CALL dbcsr_get_info(matrix_U, group=group)
432 :
433 8188 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
434 18928 : DO WHILE (dbcsr_iterator_blocks_left(iter))
435 10740 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
436 10740 : iatom = arow; CPASSERT(arow == acol)
437 10740 : CALL dbcsr_get_block_p(matrix=pao%matrix_R, row=iatom, col=iatom, block=block_R, found=found)
438 10740 : CALL dbcsr_get_block_p(matrix=matrix_U, row=iatom, col=iatom, block=block_U, found=found)
439 10740 : CPASSERT(ASSOCIATED(block_R) .AND. ASSOCIATED(block_U))
440 10740 : n = SIZE(block_U, 1)
441 :
442 : ! calculate potential V
443 32220 : ALLOCATE (vec_V(n*n))
444 646721 : vec_V(:) = 0.0_dp
445 10740 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=block_V_terms, found=found)
446 10740 : CPASSERT(ASSOCIATED(block_V_terms))
447 10740 : nterms = SIZE(block_V_terms, 2)
448 10740 : IF (nterms > 0) & ! protect against corner-case of zero pao parameters
449 71686942 : vec_V = MATMUL(block_V_terms, block_X(:, 1))
450 10740 : block_V(1:n, 1:n) => vec_V(:) ! map vector into matrix
451 :
452 : ! symmetrize
453 1430920 : IF (MAXVAL(ABS(block_V - TRANSPOSE(block_V))/MAX(1.0_dp, MAXVAL(ABS(block_V)))) > 1e-12) &
454 0 : CPABORT("block_V not symmetric")
455 1441660 : block_V = 0.5_dp*(block_V + TRANSPOSE(block_V)) ! symmetrize exactly
456 :
457 : ! regularization energy
458 : ! protect against corner-case of zero pao parameters
459 10740 : IF (PRESENT(penalty) .AND. nterms > 0) &
460 30333516 : regu_energy = regu_energy + DOT_PRODUCT(block_X(:, 1), MATMUL(block_R, block_X(:, 1)))
461 :
462 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
463 10740 : gap=gaps(iatom), evals=evals(:, iatom))
464 :
465 10740 : IF (PRESENT(matrix_G)) THEN ! TURNING POINT (if calc grad) --------------------------------
466 2132 : CPASSERT(PRESENT(matrix_M))
467 2132 : CALL dbcsr_get_block_p(matrix=matrix_M, row=iatom, col=iatom, block=block_M1, found=found)
468 :
469 : ! corner-cases: block_M1 might have been filtered out or there might be zero pao parameters
470 6396 : IF (ASSOCIATED(block_M1) .AND. SIZE(block_V_terms) > 0) THEN
471 4038 : ALLOCATE (vec_M2(n*n))
472 2019 : block_M2(1:n, 1:n) => vec_M2(:) ! map vector into matrix
473 : !TODO: this 2nd call does double work. However, *sometimes* this branch is not taken.
474 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
475 2019 : M1=block_M1, G=block_M2, gap=gaps(iatom), evals=evals(:, iatom))
476 124173 : IF (MAXVAL(ABS(block_M2 - TRANSPOSE(block_M2))) > 1e-14_dp) &
477 0 : CPABORT("matrix not symmetric")
478 :
479 : ! gradient dE/dX
480 2019 : IF (PRESENT(matrix_G)) THEN
481 2019 : CALL dbcsr_get_block_p(matrix=matrix_G, row=iatom, col=iatom, block=block_G, found=found)
482 2019 : CPASSERT(ASSOCIATED(block_G))
483 6582106 : block_G(:, 1) = MATMUL(vec_M2, block_V_terms)
484 2019 : IF (PRESENT(penalty)) &
485 10731904 : block_G = block_G + 2.0_dp*MATMUL(block_R, block_X) ! regularization gradient
486 : END IF
487 :
488 : ! forced dE/dR
489 2019 : IF (PRESENT(forces)) THEN
490 170 : ALLOCATE (M_blocks(n, n, nterms))
491 296 : DO kterm = 1, nterms
492 16806 : M_blocks(:, :, kterm) = block_M2*block_X(kterm, 1)
493 : END DO
494 34 : CALL linpot_calc_forces(pao, qs_env, iatom=iatom, M_blocks=M_blocks, forces=forces)
495 34 : DEALLOCATE (M_blocks)
496 : END IF
497 :
498 2019 : DEALLOCATE (vec_M2)
499 : END IF
500 : END IF
501 51148 : DEALLOCATE (vec_V)
502 : END DO
503 8188 : CALL dbcsr_iterator_stop(iter)
504 :
505 8188 : IF (PRESENT(penalty)) THEN
506 : ! sum penalty energies across ranks
507 7920 : CALL group%sum(penalty)
508 7920 : CALL group%sum(regu_energy)
509 7920 : penalty = penalty + regu_energy
510 : END IF
511 :
512 : ! print stuff, but not during second invocation for forces
513 8188 : IF (.NOT. PRESENT(forces)) THEN
514 : ! print eigenvalues from fock-layer
515 8154 : CALL group%sum(evals)
516 8154 : IF (pao%iw_fockev > 0) THEN
517 2000 : DO iatom = 1, natoms
518 2000 : WRITE (pao%iw_fockev, *) "PAO| atom:", iatom, " fock evals around gap:", evals(:, iatom)
519 : END DO
520 500 : CALL m_flush(pao%iw_fockev)
521 : END IF
522 : ! print homo-lumo gap encountered by fock-layer
523 8154 : CALL group%min(gaps)
524 8154 : IF (pao%iw_gap > 0) THEN
525 2000 : DO iatom = 1, natoms
526 2000 : WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
527 : END DO
528 : END IF
529 : ! one-line summaries
530 8154 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| linpot regularization energy:", regu_energy
531 37718 : IF (pao%iw > 0) WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", MINVAL(gaps), " for atom:", MINLOC(gaps)
532 : END IF
533 :
534 8188 : DEALLOCATE (gaps, evals)
535 8188 : CALL timestop(handle)
536 :
537 16376 : END SUBROUTINE pao_calc_U_linpot
538 :
539 : ! **************************************************************************************************
540 : !> \brief Internal routine, calculates terms in potential parametrization
541 : !> \param pao ...
542 : !> \param qs_env ...
543 : !> \param iatom ...
544 : !> \param V_blocks ...
545 : ! **************************************************************************************************
546 234 : SUBROUTINE linpot_calc_terms(pao, qs_env, iatom, V_blocks)
547 : TYPE(pao_env_type), POINTER :: pao
548 : TYPE(qs_environment_type), POINTER :: qs_env
549 : INTEGER, INTENT(IN) :: iatom
550 : REAL(dp), DIMENSION(:, :, :), INTENT(OUT) :: V_blocks
551 :
552 273 : SELECT CASE (pao%parameterization)
553 : CASE (pao_fock_param)
554 39 : CALL linpot_full_calc_terms(V_blocks)
555 : CASE (pao_rotinv_param)
556 195 : CALL linpot_rotinv_calc_terms(qs_env, iatom, V_blocks)
557 : CASE DEFAULT
558 234 : CPABORT("unknown parameterization")
559 : END SELECT
560 :
561 234 : END SUBROUTINE linpot_calc_terms
562 :
563 : ! **************************************************************************************************
564 : !> \brief Internal routine, calculates force contributions from potential parametrization
565 : !> \param pao ...
566 : !> \param qs_env ...
567 : !> \param iatom ...
568 : !> \param M_blocks ...
569 : !> \param forces ...
570 : ! **************************************************************************************************
571 34 : SUBROUTINE linpot_calc_forces(pao, qs_env, iatom, M_blocks, forces)
572 : TYPE(pao_env_type), POINTER :: pao
573 : TYPE(qs_environment_type), POINTER :: qs_env
574 : INTEGER, INTENT(IN) :: iatom
575 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: M_blocks
576 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: forces
577 :
578 66 : SELECT CASE (pao%parameterization)
579 : CASE (pao_fock_param)
580 : ! no force contributions
581 : CASE (pao_rotinv_param)
582 32 : CALL linpot_rotinv_calc_forces(qs_env, iatom, M_blocks, forces)
583 : CASE DEFAULT
584 34 : CPABORT("unknown parameterization")
585 : END SELECT
586 :
587 34 : END SUBROUTINE linpot_calc_forces
588 :
589 : ! **************************************************************************************************
590 : !> \brief Calculate initial guess for matrix_X
591 : !> \param pao ...
592 : !> \param qs_env ...
593 : ! **************************************************************************************************
594 34 : SUBROUTINE pao_param_initguess_linpot(pao, qs_env)
595 : TYPE(pao_env_type), POINTER :: pao
596 : TYPE(qs_environment_type), POINTER :: qs_env
597 :
598 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_initguess_linpot'
599 :
600 : INTEGER :: acol, arow, handle, i, iatom, j, k, n, &
601 : nterms
602 34 : INTEGER, DIMENSION(:), POINTER :: pri_basis_size
603 : LOGICAL :: found
604 : REAL(dp) :: w
605 34 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
606 34 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs, S_inv
607 34 : REAL(dp), DIMENSION(:), POINTER :: V_guess_vec
608 34 : REAL(dp), DIMENSION(:, :), POINTER :: block_X, V_guess, V_terms
609 : TYPE(dbcsr_iterator_type) :: iter
610 :
611 34 : CALL timeset(routineN, handle)
612 :
613 34 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis_size)
614 :
615 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,pri_basis_size) &
616 34 : !$OMP PRIVATE(iter,arow,acol,iatom,block_X,N,nterms,V_terms,found,V_guess,V_guess_vec,S,S_evecs,S_evals,S_inv,k,i,j,w)
617 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
618 : DO WHILE (dbcsr_iterator_blocks_left(iter))
619 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
620 : iatom = arow; CPASSERT(arow == acol)
621 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=V_terms, found=found)
622 : CPASSERT(ASSOCIATED(V_terms))
623 : nterms = SIZE(V_terms, 2)
624 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
625 :
626 : ! guess initial potential
627 : N = pri_basis_size(iatom)
628 : ALLOCATE (V_guess_vec(n*n))
629 : V_guess(1:n, 1:n) => V_guess_vec
630 : CALL pao_guess_initial_potential(qs_env, iatom, V_guess)
631 :
632 : ! build overlap matrix
633 : ALLOCATE (S(nterms, nterms))
634 : S(:, :) = MATMUL(TRANSPOSE(V_terms), V_terms)
635 :
636 : ! diagonalize S
637 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
638 : S_evecs(:, :) = S
639 : CALL diamat_all(S_evecs, S_evals)
640 :
641 : ! calculate Tikhonov regularized inverse
642 : ALLOCATE (S_inv(nterms, nterms))
643 : S_inv(:, :) = 0.0_dp
644 : DO k = 1, nterms
645 : w = S_evals(k)/(S_evals(k)**2 + pao%linpot_init_delta)
646 : DO i = 1, nterms
647 : DO j = 1, nterms
648 : S_inv(i, j) = S_inv(i, j) + w*S_evecs(i, k)*S_evecs(j, k)
649 : END DO
650 : END DO
651 : END DO
652 :
653 : ! perform fit
654 : block_X(:, 1) = MATMUL(MATMUL(S_inv, TRANSPOSE(V_terms)), V_guess_vec)
655 :
656 : ! clean up
657 : DEALLOCATE (V_guess_vec, S, S_evecs, S_evals, S_inv)
658 : END DO
659 : CALL dbcsr_iterator_stop(iter)
660 : !$OMP END PARALLEL
661 :
662 34 : CALL timestop(handle)
663 68 : END SUBROUTINE pao_param_initguess_linpot
664 :
665 24192 : END MODULE pao_param_linpot
|
