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