(git:f2099e5)
Loading...
Searching...
No Matches
pao_param_linpot.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
53CONTAINS
54
55! **************************************************************************************************
56!> \brief Initialize the linear potential parametrization
57!> \param pao ...
58!> \param qs_env ...
59! **************************************************************************************************
60 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 INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, row_blk_size
69 REAL(dp), DIMENSION(:, :), POINTER :: block_v_terms
70 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
75
76 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 natom=natoms)
83
84 IF (dft_control%nspins /= 1) cpabort("open shell not yet implemented")
85
86 ! figure out number of potential terms
87 ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
88 DO iatom = 1, natoms
89 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
90 CALL pao_param_count_linpot(pao, qs_env, ikind, nterms)
91 col_blk_size(iatom) = nterms
92 END DO
93
94 ! allocate matrix_V_terms
95 CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri)
96 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 col_blk_size=col_blk_size)
103 CALL dbcsr_reserve_diag_blocks(pao%matrix_V_terms)
104 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!$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 CALL pao_param_linpot_regularizer(pao)
126
127 IF (pao%precondition) THEN
128 CALL pao_param_linpot_preconditioner(pao)
129 END IF
130
131 CALL para_env%sync() ! ensure that timestop is not called too early
132
133 CALL timestop(handle)
134 END SUBROUTINE pao_param_init_linpot
135
136! **************************************************************************************************
137!> \brief Builds the regularization metric matrix_R
138!> \param pao ...
139! **************************************************************************************************
140 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 INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
148 LOGICAL :: found
149 REAL(dp) :: v, w
150 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
151 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs
152 REAL(dp), DIMENSION(:, :), POINTER :: block_r, v_terms
153 TYPE(dbcsr_iterator_type) :: iter
154
155 CALL timeset(routinen, handle)
156
157 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot regularizer"
158
159 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 name="PAO matrix_R")
168 CALL dbcsr_reserve_diag_blocks(pao%matrix_R)
169
170 ! fill matrix_R
171!$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
172!$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 CALL timestop(handle)
209 END SUBROUTINE pao_param_linpot_regularizer
210
211! **************************************************************************************************
212!> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
213!> \param pao ...
214! **************************************************************************************************
215 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 INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
223 LOGICAL :: found
224 REAL(dp) :: eval_capped
225 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
226 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs
227 REAL(dp), DIMENSION(:, :), POINTER :: block_precon, block_precon_inv, &
228 block_v_terms
229 TYPE(dbcsr_iterator_type) :: iter
230
231 CALL timeset(routinen, handle)
232
233 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot preconditioner"
234
235 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 name="PAO matrix_precon")
243 CALL dbcsr_reserve_diag_blocks(pao%matrix_precon)
244
245 CALL dbcsr_create(pao%matrix_precon_inv, template=pao%matrix_precon, name="PAO matrix_precon_inv")
246 CALL dbcsr_reserve_diag_blocks(pao%matrix_precon_inv)
247
248!$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
249!$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 CALL timestop(handle)
289 END SUBROUTINE pao_param_linpot_preconditioner
290
291! **************************************************************************************************
292!> \brief Finalize the linear potential parametrization
293!> \param pao ...
294! **************************************************************************************************
295 SUBROUTINE pao_param_finalize_linpot(pao)
296 TYPE(pao_env_type), POINTER :: pao
297
298 CALL dbcsr_release(pao%matrix_V_terms)
299 CALL dbcsr_release(pao%matrix_R)
300
301 IF (pao%precondition) THEN
302 CALL dbcsr_release(pao%matrix_precon)
303 CALL dbcsr_release(pao%matrix_precon_inv)
304 END IF
305
306 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 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
324
325 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 pao_basis_size=pao_basis_size)
330
331 IF (pao_basis_size == basis_set%nsgf) THEN
332 nparams = 0 ! pao disabled for iatom
333
334 ELSE
335 SELECT CASE (pao%parameterization)
336 CASE (pao_fock_param)
337 CALL linpot_full_count_terms(qs_env, ikind, nterms=nparams)
338 CASE (pao_rotinv_param)
339 CALL linpot_rotinv_count_terms(qs_env, ikind, nterms=nparams)
340 CASE DEFAULT
341 cpabort("unknown parameterization")
342 END SELECT
343 END IF
344
345 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
368 TYPE(dbcsr_type) :: matrix_m, matrix_u
369
370 CALL timeset(routinen, handle)
371 CALL get_qs_env(qs_env, matrix_s=matrix_s)
372 CALL dbcsr_create(matrix_u, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
373 CALL dbcsr_reserve_diag_blocks(matrix_u)
374
375 !TODO: move this condition into pao_calc_U, use matrix_N as template
376 IF (gradient) THEN
377 CALL pao_calc_grad_lnv_wrt_u(qs_env, ls_scf_env, matrix_m)
378 CALL pao_calc_u_linpot(pao, qs_env, matrix_u, matrix_m, pao%matrix_G, penalty, forces)
379 CALL dbcsr_release(matrix_m)
380 ELSE
381 CALL pao_calc_u_linpot(pao, qs_env, matrix_u, penalty=penalty)
382 END IF
383
384 CALL pao_calc_ab_from_u(pao, qs_env, ls_scf_env, matrix_u)
385 CALL dbcsr_release(matrix_u)
386 CALL timestop(handle)
387 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 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 REAL(kind=dp), PARAMETER :: eps_block_v_sym = 1.0e-12_dp
409
410 INTEGER :: acol, arow, handle, iatom, kterm, n, &
411 natoms, nterms
412 LOGICAL :: found
413 REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
414 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: evals
415 REAL(dp), DIMENSION(:), POINTER :: vec_m2, vec_v
416 REAL(dp), DIMENSION(:, :), POINTER :: block_g, block_m1, block_m2, block_r, &
417 block_u, block_v, block_v_terms, &
418 block_x
419 REAL(dp), DIMENSION(:, :, :), POINTER :: m_blocks
420 REAL(kind=dp) :: regu_energy
421 TYPE(dbcsr_iterator_type) :: iter
422 TYPE(mp_comm_type) :: group
423
424 CALL timeset(routinen, handle)
425
426 cpassert(PRESENT(matrix_g) .EQV. PRESENT(matrix_m))
427
428 CALL get_qs_env(qs_env, natom=natoms)
429 ALLOCATE (gaps(natoms), evals(10, natoms)) ! printing 10 eigenvalues seems reasonable
430 evals(:, :) = 0.0_dp
431 gaps(:) = huge(1.0_dp)
432 regu_energy = 0.0_dp
433 CALL dbcsr_get_info(matrix_u, group=group)
434
435 CALL dbcsr_iterator_start(iter, pao%matrix_X)
436 DO WHILE (dbcsr_iterator_blocks_left(iter))
437 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
438 iatom = arow; cpassert(arow == acol)
439 CALL dbcsr_get_block_p(matrix=pao%matrix_R, row=iatom, col=iatom, block=block_r, found=found)
440 CALL dbcsr_get_block_p(matrix=matrix_u, row=iatom, col=iatom, block=block_u, found=found)
441 cpassert(ASSOCIATED(block_r) .AND. ASSOCIATED(block_u))
442 n = SIZE(block_u, 1)
443
444 ! calculate potential V
445 ALLOCATE (vec_v(n*n))
446 vec_v(:) = 0.0_dp
447 CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=block_v_terms, found=found)
448 cpassert(ASSOCIATED(block_v_terms))
449 nterms = SIZE(block_v_terms, 2)
450 IF (nterms > 0) THEN
451 ! protect against corner-case of zero pao parameters
452 vec_v = matmul(block_v_terms, block_x(:, 1))
453 END IF
454 block_v(1:n, 1:n) => vec_v(:) ! map vector into matrix
455
456 ! symmetrize
457 IF (maxval(abs(block_v - transpose(block_v))/max(1.0_dp, maxval(abs(block_v)))) > eps_block_v_sym) THEN
458 cpabort("block_V not symmetric")
459 END IF
460 block_v = 0.5_dp*(block_v + transpose(block_v)) ! symmetrize exactly
461
462 ! regularization energy
463 ! protect against corner-case of zero pao parameters
464 IF (PRESENT(penalty) .AND. nterms > 0) THEN
465 regu_energy = regu_energy + dot_product(block_x(:, 1), matmul(block_r, block_x(:, 1)))
466 END IF
467
468 CALL pao_calc_u_block_fock(pao, iatom=iatom, penalty=penalty, v=block_v, u=block_u, &
469 gap=gaps(iatom), evals=evals(:, iatom))
470
471 IF (PRESENT(matrix_g)) THEN ! TURNING POINT (if calc grad) --------------------------------
472 cpassert(PRESENT(matrix_m))
473 CALL dbcsr_get_block_p(matrix=matrix_m, row=iatom, col=iatom, block=block_m1, found=found)
474
475 ! corner-cases: block_M1 might have been filtered out or there might be zero pao parameters
476 IF (ASSOCIATED(block_m1) .AND. SIZE(block_v_terms) > 0) THEN
477 ALLOCATE (vec_m2(n*n))
478 block_m2(1:n, 1:n) => vec_m2(:) ! map vector into matrix
479 !TODO: this 2nd call does double work. However, *sometimes* this branch is not taken.
480 CALL pao_calc_u_block_fock(pao, iatom=iatom, penalty=penalty, v=block_v, u=block_u, &
481 m1=block_m1, g=block_m2, gap=gaps(iatom), evals=evals(:, iatom))
482 IF (maxval(abs(block_m2 - transpose(block_m2))) > 1e-14_dp) THEN
483 cpabort("matrix not symmetric")
484 END IF
485
486 ! gradient dE/dX
487 IF (PRESENT(matrix_g)) THEN
488 CALL dbcsr_get_block_p(matrix=matrix_g, row=iatom, col=iatom, block=block_g, found=found)
489 cpassert(ASSOCIATED(block_g))
490 block_g(:, 1) = matmul(vec_m2, block_v_terms)
491 IF (PRESENT(penalty)) THEN
492 block_g = block_g + 2.0_dp*matmul(block_r, block_x)
493 END IF ! regularization gradient
494 END IF
495
496 ! forced dE/dR
497 IF (PRESENT(forces)) THEN
498 ALLOCATE (m_blocks(n, n, nterms))
499 DO kterm = 1, nterms
500 m_blocks(:, :, kterm) = block_m2*block_x(kterm, 1)
501 END DO
502 CALL linpot_calc_forces(pao, qs_env, iatom=iatom, m_blocks=m_blocks, forces=forces)
503 DEALLOCATE (m_blocks)
504 END IF
505
506 DEALLOCATE (vec_m2)
507 END IF
508 END IF
509 DEALLOCATE (vec_v)
510 END DO
511 CALL dbcsr_iterator_stop(iter)
512
513 IF (PRESENT(penalty)) THEN
514 ! sum penalty energies across ranks
515 CALL group%sum(penalty)
516 CALL group%sum(regu_energy)
517 penalty = penalty + regu_energy
518 END IF
519
520 ! print stuff, but not during second invocation for forces
521 IF (.NOT. PRESENT(forces)) THEN
522 ! print eigenvalues from fock-layer
523 CALL group%sum(evals)
524 IF (pao%iw_fockev > 0) THEN
525 DO iatom = 1, natoms
526 WRITE (pao%iw_fockev, *) "PAO| atom:", iatom, " fock evals around gap:", evals(:, iatom)
527 END DO
528 CALL m_flush(pao%iw_fockev)
529 END IF
530 ! print homo-lumo gap encountered by fock-layer
531 CALL group%min(gaps)
532 IF (pao%iw_gap > 0) THEN
533 DO iatom = 1, natoms
534 WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
535 END DO
536 END IF
537 ! one-line summaries
538 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| linpot regularization energy:", regu_energy
539 IF (pao%iw > 0) WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", minval(gaps), " for atom:", minloc(gaps)
540 END IF
541
542 DEALLOCATE (gaps, evals)
543 CALL timestop(handle)
544
545 END SUBROUTINE pao_calc_u_linpot
546
547! **************************************************************************************************
548!> \brief Internal routine, calculates terms in potential parametrization
549!> \param pao ...
550!> \param qs_env ...
551!> \param iatom ...
552!> \param V_blocks ...
553! **************************************************************************************************
554 SUBROUTINE linpot_calc_terms(pao, qs_env, iatom, V_blocks)
555 TYPE(pao_env_type), POINTER :: pao
556 TYPE(qs_environment_type), POINTER :: qs_env
557 INTEGER, INTENT(IN) :: iatom
558 REAL(dp), DIMENSION(:, :, :), INTENT(OUT) :: v_blocks
559
560 SELECT CASE (pao%parameterization)
561 CASE (pao_fock_param)
562 CALL linpot_full_calc_terms(v_blocks)
563 CASE (pao_rotinv_param)
564 CALL linpot_rotinv_calc_terms(qs_env, iatom, v_blocks)
565 CASE DEFAULT
566 cpabort("unknown parameterization")
567 END SELECT
568
569 END SUBROUTINE linpot_calc_terms
570
571! **************************************************************************************************
572!> \brief Internal routine, calculates force contributions from potential parametrization
573!> \param pao ...
574!> \param qs_env ...
575!> \param iatom ...
576!> \param M_blocks ...
577!> \param forces ...
578! **************************************************************************************************
579 SUBROUTINE linpot_calc_forces(pao, qs_env, iatom, M_blocks, forces)
580 TYPE(pao_env_type), POINTER :: pao
581 TYPE(qs_environment_type), POINTER :: qs_env
582 INTEGER, INTENT(IN) :: iatom
583 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: m_blocks
584 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: forces
585
586 SELECT CASE (pao%parameterization)
587 CASE (pao_fock_param)
588 ! no force contributions
589 CASE (pao_rotinv_param)
590 CALL linpot_rotinv_calc_forces(qs_env, iatom, m_blocks, forces)
591 CASE DEFAULT
592 cpabort("unknown parameterization")
593 END SELECT
594
595 END SUBROUTINE linpot_calc_forces
596
597! **************************************************************************************************
598!> \brief Calculate initial guess for matrix_X
599!> \param pao ...
600!> \param qs_env ...
601! **************************************************************************************************
602 SUBROUTINE pao_param_initguess_linpot(pao, qs_env)
603 TYPE(pao_env_type), POINTER :: pao
604 TYPE(qs_environment_type), POINTER :: qs_env
605
606 CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_initguess_linpot'
607
608 INTEGER :: acol, arow, handle, i, iatom, j, k, n, &
609 nterms
610 INTEGER, DIMENSION(:), POINTER :: pri_basis_size
611 LOGICAL :: found
612 REAL(dp) :: w
613 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
614 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs, s_inv
615 REAL(dp), DIMENSION(:), POINTER :: v_guess_vec
616 REAL(dp), DIMENSION(:, :), POINTER :: block_x, v_guess, v_terms
617 TYPE(dbcsr_iterator_type) :: iter
618
619 CALL timeset(routinen, handle)
620
621 CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis_size)
622
623!$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,pri_basis_size) &
624!$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)
625 CALL dbcsr_iterator_start(iter, pao%matrix_X)
626 DO WHILE (dbcsr_iterator_blocks_left(iter))
627 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
628 iatom = arow; cpassert(arow == acol)
629 CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=v_terms, found=found)
630 cpassert(ASSOCIATED(v_terms))
631 nterms = SIZE(v_terms, 2)
632 IF (nterms == 0) cycle ! protect against corner-case of zero pao parameters
633
634 ! guess initial potential
635 n = pri_basis_size(iatom)
636 ALLOCATE (v_guess_vec(n*n))
637 v_guess(1:n, 1:n) => v_guess_vec
638 CALL pao_guess_initial_potential(qs_env, iatom, v_guess)
639
640 ! build overlap matrix
641 ALLOCATE (s(nterms, nterms))
642 s(:, :) = matmul(transpose(v_terms), v_terms)
643
644 ! diagonalize S
645 ALLOCATE (s_evals(nterms), s_evecs(nterms, nterms))
646 s_evecs(:, :) = s
647 CALL diamat_all(s_evecs, s_evals)
648
649 ! calculate Tikhonov regularized inverse
650 ALLOCATE (s_inv(nterms, nterms))
651 s_inv(:, :) = 0.0_dp
652 DO k = 1, nterms
653 w = s_evals(k)/(s_evals(k)**2 + pao%linpot_init_delta)
654 DO i = 1, nterms
655 DO j = 1, nterms
656 s_inv(i, j) = s_inv(i, j) + w*s_evecs(i, k)*s_evecs(j, k)
657 END DO
658 END DO
659 END DO
660
661 ! perform fit
662 block_x(:, 1) = matmul(matmul(s_inv, transpose(v_terms)), v_guess_vec)
663
664 ! clean up
665 DEALLOCATE (v_guess_vec, s, s_evecs, s_evals, s_inv)
666 END DO
667 CALL dbcsr_iterator_stop(iter)
668!$OMP END PARALLEL
669
670 CALL timestop(handle)
671 END SUBROUTINE pao_param_initguess_linpot
672
673END MODULE pao_param_linpot
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_reserve_diag_blocks
subroutine, public dbcsr_reserve_diag_blocks(matrix)
Reserves all diagonal blocks.
Definition cp_dbcsr_contrib.F:255
dm_ls_scf_types
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Definition dm_ls_scf_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:381
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
pao_input
Definition pao_input.F:8
pao_input::pao_fock_param
integer, parameter, public pao_fock_param
Definition pao_input.F:45
pao_input::pao_rotinv_param
integer, parameter, public pao_rotinv_param
Definition pao_input.F:45
pao_linpot_full
Full parametrization of Fock matrix, ie. the identity parametrization.
Definition pao_linpot_full.F:12
pao_linpot_full::linpot_full_calc_terms
subroutine, public linpot_full_calc_terms(v_blocks)
Builds potential terms.
Definition pao_linpot_full.F:58
pao_linpot_full::linpot_full_count_terms
subroutine, public linpot_full_count_terms(qs_env, ikind, nterms)
Count number of terms for given atomic kind.
Definition pao_linpot_full.F:38
pao_linpot_rotinv
Rotationally invariant parametrization of Fock matrix.
Definition pao_linpot_rotinv.F:12
pao_linpot_rotinv::linpot_rotinv_calc_forces
subroutine, public linpot_rotinv_calc_forces(qs_env, iatom, m_blocks, forces)
Calculate force contribution from rotinv parametrization.
Definition pao_linpot_rotinv.F:336
pao_linpot_rotinv::linpot_rotinv_calc_terms
subroutine, public linpot_rotinv_calc_terms(qs_env, iatom, v_blocks)
Calculate all potential terms of the rotinv parametrization.
Definition pao_linpot_rotinv.F:123
pao_linpot_rotinv::linpot_rotinv_count_terms
subroutine, public linpot_rotinv_count_terms(qs_env, ikind, nterms)
Count number of terms for given atomic kind.
Definition pao_linpot_rotinv.F:48
pao_param_fock
Common framework for using eigenvectors of a Fock matrix as PAO basis.
Definition pao_param_fock.F:12
pao_param_fock::pao_calc_u_block_fock
subroutine, public pao_calc_u_block_fock(pao, iatom, v, u, penalty, gap, evals, m1, g)
Calculate new matrix U and optinally its gradient G.
Definition pao_param_fock.F:43
pao_param_linpot
Common framework for a linear parametrization of the potential.
Definition pao_param_linpot.F:12
pao_param_linpot::pao_param_finalize_linpot
subroutine, public pao_param_finalize_linpot(pao)
Finalize the linear potential parametrization.
Definition pao_param_linpot.F:296
pao_param_linpot::pao_param_init_linpot
subroutine, public pao_param_init_linpot(pao, qs_env)
Initialize the linear potential parametrization.
Definition pao_param_linpot.F:61
pao_param_linpot::pao_calc_ab_linpot
subroutine, public pao_calc_ab_linpot(pao, qs_env, ls_scf_env, gradient, penalty, forces)
Takes current matrix_X and calculates the matrices A and B.
Definition pao_param_linpot.F:357
pao_param_linpot::pao_param_count_linpot
subroutine, public pao_param_count_linpot(pao, qs_env, ikind, nparams)
Returns the number of potential terms for given atomic kind.
Definition pao_param_linpot.F:316
pao_param_linpot::pao_param_initguess_linpot
subroutine, public pao_param_initguess_linpot(pao, qs_env)
Calculate initial guess for matrix_X.
Definition pao_param_linpot.F:603
pao_param_methods
Common routines for PAO parametrizations.
Definition pao_param_methods.F:12
pao_param_methods::pao_calc_grad_lnv_wrt_u
subroutine, public pao_calc_grad_lnv_wrt_u(qs_env, ls_scf_env, matrix_m_diag)
Helper routine, calculates partial derivative dE/dU.
Definition pao_param_methods.F:50
pao_param_methods::pao_calc_ab_from_u
subroutine, public pao_calc_ab_from_u(pao, qs_env, ls_scf_env, matrix_u_diag)
Takes current matrix_X and calculates the matrices A and B.
Definition pao_param_methods.F:115
pao_potentials
Factory routines for potentials used e.g. by pao_param_exp and pao_ml.
Definition pao_potentials.F:12
pao_potentials::pao_guess_initial_potential
subroutine, public pao_guess_initial_potential(qs_env, iatom, block_v)
Makes an educated guess for the initial potential based on positions of neighboring atoms.
Definition pao_potentials.F:49
pao_types
Types used by the PAO machinery.
Definition pao_types.F:12
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::dft_control_type
Definition cp_control_types.F:780
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
dm_ls_scf_types::ls_scf_env_type
Definition dm_ls_scf_types.F:88
message_passing::mp_comm_type
Definition message_passing.F:179
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
pao_types::pao_env_type
Definition pao_types.F:148
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190