(git:ed6f26b)
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-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! **************************************************************************************************
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) &
128 CALL pao_param_linpot_preconditioner(pao)
129
130 CALL para_env%sync() ! ensure that timestop is not called too early
131
132 CALL timestop(handle)
133 END SUBROUTINE pao_param_init_linpot
134
135! **************************************************************************************************
136!> \brief Builds the regularization metric matrix_R
137!> \param pao ...
138! **************************************************************************************************
139 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 INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
147 LOGICAL :: found
148 REAL(dp) :: v, w
149 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
150 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs
151 REAL(dp), DIMENSION(:, :), POINTER :: block_r, v_terms
152 TYPE(dbcsr_iterator_type) :: iter
153
154 CALL timeset(routinen, handle)
155
156 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot regularizer"
157
158 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 name="PAO matrix_R")
167 CALL dbcsr_reserve_diag_blocks(pao%matrix_R)
168
169 ! fill matrix_R
170!$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
171!$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 CALL timestop(handle)
208 END SUBROUTINE pao_param_linpot_regularizer
209
210! **************************************************************************************************
211!> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
212!> \param pao ...
213! **************************************************************************************************
214 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 INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
222 LOGICAL :: found
223 REAL(dp) :: eval_capped
224 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
225 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs
226 REAL(dp), DIMENSION(:, :), POINTER :: block_precon, block_precon_inv, &
227 block_v_terms
228 TYPE(dbcsr_iterator_type) :: iter
229
230 CALL timeset(routinen, handle)
231
232 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot preconditioner"
233
234 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 name="PAO matrix_precon")
242 CALL dbcsr_reserve_diag_blocks(pao%matrix_precon)
243
244 CALL dbcsr_create(pao%matrix_precon_inv, template=pao%matrix_precon, name="PAO matrix_precon_inv")
245 CALL dbcsr_reserve_diag_blocks(pao%matrix_precon_inv)
246
247!$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
248!$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 CALL timestop(handle)
288 END SUBROUTINE pao_param_linpot_preconditioner
289
290! **************************************************************************************************
291!> \brief Finalize the linear potential parametrization
292!> \param pao ...
293! **************************************************************************************************
294 SUBROUTINE pao_param_finalize_linpot(pao)
295 TYPE(pao_env_type), POINTER :: pao
296
297 CALL dbcsr_release(pao%matrix_V_terms)
298 CALL dbcsr_release(pao%matrix_R)
299
300 IF (pao%precondition) THEN
301 CALL dbcsr_release(pao%matrix_precon)
302 CALL dbcsr_release(pao%matrix_precon_inv)
303 END IF
304
305 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 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
323
324 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 pao_basis_size=pao_basis_size)
329
330 IF (pao_basis_size == basis_set%nsgf) THEN
331 nparams = 0 ! pao disabled for iatom
332
333 ELSE
334 SELECT CASE (pao%parameterization)
335 CASE (pao_fock_param)
336 CALL linpot_full_count_terms(qs_env, ikind, nterms=nparams)
337 CASE (pao_rotinv_param)
338 CALL linpot_rotinv_count_terms(qs_env, ikind, nterms=nparams)
339 CASE DEFAULT
340 cpabort("unknown parameterization")
341 END SELECT
342 END IF
343
344 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
367 TYPE(dbcsr_type) :: matrix_m, matrix_u
368
369 CALL timeset(routinen, handle)
370 CALL get_qs_env(qs_env, matrix_s=matrix_s)
371 CALL dbcsr_create(matrix_u, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
372 CALL dbcsr_reserve_diag_blocks(matrix_u)
373
374 !TODO: move this condition into pao_calc_U, use matrix_N as template
375 IF (gradient) THEN
376 CALL pao_calc_grad_lnv_wrt_u(qs_env, ls_scf_env, matrix_m)
377 CALL pao_calc_u_linpot(pao, qs_env, matrix_u, matrix_m, pao%matrix_G, penalty, forces)
378 CALL dbcsr_release(matrix_m)
379 ELSE
380 CALL pao_calc_u_linpot(pao, qs_env, matrix_u, penalty=penalty)
381 END IF
382
383 CALL pao_calc_ab_from_u(pao, qs_env, ls_scf_env, matrix_u)
384 CALL dbcsr_release(matrix_u)
385 CALL timestop(handle)
386 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 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 REAL(dp), DIMENSION(:), POINTER :: vec_m2, vec_v
414 REAL(dp), DIMENSION(:, :), POINTER :: block_g, block_m1, block_m2, block_r, &
415 block_u, block_v, block_v_terms, &
416 block_x
417 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 CALL timeset(routinen, handle)
423
424 cpassert(PRESENT(matrix_g) .EQV. PRESENT(matrix_m))
425
426 CALL get_qs_env(qs_env, natom=natoms)
427 ALLOCATE (gaps(natoms), evals(10, natoms)) ! printing 10 eigenvalues seems reasonable
428 evals(:, :) = 0.0_dp
429 gaps(:) = huge(1.0_dp)
430 regu_energy = 0.0_dp
431 CALL dbcsr_get_info(matrix_u, group=group)
432
433 CALL dbcsr_iterator_start(iter, pao%matrix_X)
434 DO WHILE (dbcsr_iterator_blocks_left(iter))
435 CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
436 iatom = arow; cpassert(arow == acol)
437 CALL dbcsr_get_block_p(matrix=pao%matrix_R, row=iatom, col=iatom, block=block_r, found=found)
438 CALL dbcsr_get_block_p(matrix=matrix_u, row=iatom, col=iatom, block=block_u, found=found)
439 cpassert(ASSOCIATED(block_r) .AND. ASSOCIATED(block_u))
440 n = SIZE(block_u, 1)
441
442 ! calculate potential V
443 ALLOCATE (vec_v(n*n))
444 vec_v(:) = 0.0_dp
445 CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=block_v_terms, found=found)
446 cpassert(ASSOCIATED(block_v_terms))
447 nterms = SIZE(block_v_terms, 2)
448 IF (nterms > 0) & ! protect against corner-case of zero pao parameters
449 vec_v = matmul(block_v_terms, block_x(:, 1))
450 block_v(1:n, 1:n) => vec_v(:) ! map vector into matrix
451
452 ! symmetrize
453 IF (maxval(abs(block_v - transpose(block_v))/max(1.0_dp, maxval(abs(block_v)))) > 1e-12) &
454 cpabort("block_V not symmetric")
455 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 IF (PRESENT(penalty) .AND. nterms > 0) &
460 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 gap=gaps(iatom), evals=evals(:, iatom))
464
465 IF (PRESENT(matrix_g)) THEN ! TURNING POINT (if calc grad) --------------------------------
466 cpassert(PRESENT(matrix_m))
467 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 IF (ASSOCIATED(block_m1) .AND. SIZE(block_v_terms) > 0) THEN
471 ALLOCATE (vec_m2(n*n))
472 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 m1=block_m1, g=block_m2, gap=gaps(iatom), evals=evals(:, iatom))
476 IF (maxval(abs(block_m2 - transpose(block_m2))) > 1e-14_dp) &
477 cpabort("matrix not symmetric")
478
479 ! gradient dE/dX
480 IF (PRESENT(matrix_g)) THEN
481 CALL dbcsr_get_block_p(matrix=matrix_g, row=iatom, col=iatom, block=block_g, found=found)
482 cpassert(ASSOCIATED(block_g))
483 block_g(:, 1) = matmul(vec_m2, block_v_terms)
484 IF (PRESENT(penalty)) &
485 block_g = block_g + 2.0_dp*matmul(block_r, block_x) ! regularization gradient
486 END IF
487
488 ! forced dE/dR
489 IF (PRESENT(forces)) THEN
490 ALLOCATE (m_blocks(n, n, nterms))
491 DO kterm = 1, nterms
492 m_blocks(:, :, kterm) = block_m2*block_x(kterm, 1)
493 END DO
494 CALL linpot_calc_forces(pao, qs_env, iatom=iatom, m_blocks=m_blocks, forces=forces)
495 DEALLOCATE (m_blocks)
496 END IF
497
498 DEALLOCATE (vec_m2)
499 END IF
500 END IF
501 DEALLOCATE (vec_v)
502 END DO
503 CALL dbcsr_iterator_stop(iter)
504
505 IF (PRESENT(penalty)) THEN
506 ! sum penalty energies across ranks
507 CALL group%sum(penalty)
508 CALL group%sum(regu_energy)
509 penalty = penalty + regu_energy
510 END IF
511
512 ! print stuff, but not during second invocation for forces
513 IF (.NOT. PRESENT(forces)) THEN
514 ! print eigenvalues from fock-layer
515 CALL group%sum(evals)
516 IF (pao%iw_fockev > 0) THEN
517 DO iatom = 1, natoms
518 WRITE (pao%iw_fockev, *) "PAO| atom:", iatom, " fock evals around gap:", evals(:, iatom)
519 END DO
520 CALL m_flush(pao%iw_fockev)
521 END IF
522 ! print homo-lumo gap encountered by fock-layer
523 CALL group%min(gaps)
524 IF (pao%iw_gap > 0) THEN
525 DO iatom = 1, natoms
526 WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
527 END DO
528 END IF
529 ! one-line summaries
530 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| linpot regularization energy:", regu_energy
531 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 DEALLOCATE (gaps, evals)
535 CALL timestop(handle)
536
537 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 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 SELECT CASE (pao%parameterization)
553 CASE (pao_fock_param)
554 CALL linpot_full_calc_terms(v_blocks)
555 CASE (pao_rotinv_param)
556 CALL linpot_rotinv_calc_terms(qs_env, iatom, v_blocks)
557 CASE DEFAULT
558 cpabort("unknown parameterization")
559 END SELECT
560
561 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 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 SELECT CASE (pao%parameterization)
579 CASE (pao_fock_param)
580 ! no force contributions
581 CASE (pao_rotinv_param)
582 CALL linpot_rotinv_calc_forces(qs_env, iatom, m_blocks, forces)
583 CASE DEFAULT
584 cpabort("unknown parameterization")
585 END SELECT
586
587 END SUBROUTINE linpot_calc_forces
588
589! **************************************************************************************************
590!> \brief Calculate initial guess for matrix_X
591!> \param pao ...
592!> \param qs_env ...
593! **************************************************************************************************
594 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 INTEGER, DIMENSION(:), POINTER :: pri_basis_size
603 LOGICAL :: found
604 REAL(dp) :: w
605 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
606 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs, s_inv
607 REAL(dp), DIMENSION(:), POINTER :: v_guess_vec
608 REAL(dp), DIMENSION(:, :), POINTER :: block_x, v_guess, v_terms
609 TYPE(dbcsr_iterator_type) :: iter
610
611 CALL timeset(routinen, handle)
612
613 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!$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 CALL timestop(handle)
663 END SUBROUTINE pao_param_initguess_linpot
664
665END MODULE pao_param_linpot
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
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_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)
...
Definition cp_dbcsr_api.F:956
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:930
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1027
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:679
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:794
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:989
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
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:130
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:373
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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:333
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:121
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:295
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:356
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:315
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:595
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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, init_u_ramping_each_scf, reltmat, ghost, 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:454
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:186
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
dm_ls_scf_types::ls_scf_env_type
Definition dm_ls_scf_types.F:88
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pao_types::pao_env_type
Definition pao_types.F:147
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171