(git:374b731)
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-2024 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 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_reserve_diag_blocks, dbcsr_type
20 USE dm_ls_scf_types, ONLY: ls_scf_env_type
21 USE kinds, ONLY: dp
22 USE machine, ONLY: m_flush
23 USE mathlib, ONLY: diamat_all
24 USE message_passing, ONLY: mp_comm_type,&
25 mp_para_env_type
26 USE pao_input, ONLY: pao_fock_param,&
27 pao_rotinv_param
28 USE pao_linpot_full, ONLY: linpot_full_calc_terms,&
29 linpot_full_count_terms
30 USE pao_linpot_rotinv, ONLY: linpot_rotinv_calc_forces,&
31 linpot_rotinv_calc_terms,&
32 linpot_rotinv_count_terms
33 USE pao_param_fock, ONLY: pao_calc_u_block_fock
34 USE pao_param_methods, ONLY: pao_calc_ab_from_u,&
35 pao_calc_grad_lnv_wrt_u
36 USE pao_potentials, ONLY: pao_guess_initial_potential
37 USE pao_types, ONLY: pao_env_type
38 USE particle_types, ONLY: particle_type
39 USE qs_environment_types, ONLY: get_qs_env,&
40 qs_environment_type
41 USE qs_kind_types, ONLY: get_qs_kind,&
42 qs_kind_type
43#include "./base/base_uses.f90"
44
45 IMPLICIT NONE
46
47 PRIVATE
48
49 PUBLIC :: pao_param_init_linpot, pao_param_finalize_linpot, pao_calc_ab_linpot
50 PUBLIC :: pao_param_count_linpot, pao_param_initguess_linpot
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief Initialize the linear potential parametrization
56!> \param pao ...
57!> \param qs_env ...
58! **************************************************************************************************
59 SUBROUTINE pao_param_init_linpot(pao, qs_env)
60 TYPE(pao_env_type), POINTER :: pao
61 TYPE(qs_environment_type), POINTER :: qs_env
62
63 CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_init_linpot'
64
65 INTEGER :: acol, arow, handle, iatom, ikind, n, &
66 natoms, nterms
67 INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, row_blk_size
68 REAL(dp), DIMENSION(:, :), POINTER :: block_v_terms
69 REAL(dp), DIMENSION(:, :, :), POINTER :: v_blocks
70 TYPE(dbcsr_iterator_type) :: iter
71 TYPE(dft_control_type), POINTER :: dft_control
72 TYPE(mp_para_env_type), POINTER :: para_env
73 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
74
75 CALL timeset(routinen, handle)
76
77 CALL get_qs_env(qs_env, &
78 para_env=para_env, &
79 dft_control=dft_control, &
80 particle_set=particle_set, &
81 natom=natoms)
82
83 IF (dft_control%nspins /= 1) cpabort("open shell not yet implemented")
84
85 ! figure out number of potential terms
86 ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
87 DO iatom = 1, natoms
88 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
89 CALL pao_param_count_linpot(pao, qs_env, ikind, nterms)
90 col_blk_size(iatom) = nterms
91 END DO
92
93 ! allocate matrix_V_terms
94 CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri)
95 row_blk_size = blk_sizes_pri**2
96 CALL dbcsr_create(pao%matrix_V_terms, &
97 name="PAO matrix_V_terms", &
98 dist=pao%diag_distribution, &
99 matrix_type="N", &
100 row_blk_size=row_blk_size, &
101 col_blk_size=col_blk_size)
102 CALL dbcsr_reserve_diag_blocks(pao%matrix_V_terms)
103 DEALLOCATE (row_blk_size, col_blk_size)
104
105 ! calculate, normalize, and store potential terms as rows of block_V_terms
106!$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,blk_sizes_pri) &
107!$OMP PRIVATE(iter,arow,acol,iatom,N,nterms,block_V_terms,V_blocks)
108 CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
109 DO WHILE (dbcsr_iterator_blocks_left(iter))
110 CALL dbcsr_iterator_next_block(iter, arow, acol, block_v_terms)
111 iatom = arow; cpassert(arow == acol)
112 nterms = SIZE(block_v_terms, 2)
113 IF (nterms == 0) cycle ! protect against corner-case of zero pao parameters
114 n = blk_sizes_pri(iatom)
115 cpassert(n*n == SIZE(block_v_terms, 1))
116 ALLOCATE (v_blocks(n, n, nterms))
117 CALL linpot_calc_terms(pao, qs_env, iatom, v_blocks)
118 block_v_terms = reshape(v_blocks, (/n*n, nterms/)) ! convert matrices into vectors
119 DEALLOCATE (v_blocks)
120 END DO
121 CALL dbcsr_iterator_stop(iter)
122!$OMP END PARALLEL
123
124 CALL pao_param_linpot_regularizer(pao)
125
126 IF (pao%precondition) &
127 CALL pao_param_linpot_preconditioner(pao)
128
129 CALL para_env%sync() ! ensure that timestop is not called too early
130
131 CALL timestop(handle)
132 END SUBROUTINE pao_param_init_linpot
133
134! **************************************************************************************************
135!> \brief Builds the regularization metric matrix_R
136!> \param pao ...
137! **************************************************************************************************
138 SUBROUTINE pao_param_linpot_regularizer(pao)
139 TYPE(pao_env_type), POINTER :: pao
140
141 CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_linpot_regularizer'
142
143 INTEGER :: acol, arow, handle, i, iatom, j, k, &
144 nterms
145 INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
146 LOGICAL :: found
147 REAL(dp) :: v, w
148 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
149 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs
150 REAL(dp), DIMENSION(:, :), POINTER :: block_r, v_terms
151 TYPE(dbcsr_iterator_type) :: iter
152
153 CALL timeset(routinen, handle)
154
155 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot regularizer"
156
157 CALL dbcsr_get_info(pao%matrix_V_terms, col_blk_size=blk_sizes_nterms)
158
159 ! build regularization metric
160 CALL dbcsr_create(pao%matrix_R, &
161 template=pao%matrix_V_terms, &
162 matrix_type="N", &
163 row_blk_size=blk_sizes_nterms, &
164 col_blk_size=blk_sizes_nterms, &
165 name="PAO matrix_R")
166 CALL dbcsr_reserve_diag_blocks(pao%matrix_R)
167
168 ! fill matrix_R
169!$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
170!$OMP PRIVATE(iter,arow,acol,iatom,block_R,V_terms,found,nterms,S,S_evecs,S_evals,k,i,j,v,w)
171 CALL dbcsr_iterator_start(iter, pao%matrix_R)
172 DO WHILE (dbcsr_iterator_blocks_left(iter))
173 CALL dbcsr_iterator_next_block(iter, arow, acol, block_r)
174 iatom = arow; cpassert(arow == acol)
175 CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=v_terms, found=found)
176 cpassert(ASSOCIATED(v_terms))
177 nterms = SIZE(v_terms, 2)
178 IF (nterms == 0) cycle ! protect against corner-case of zero pao parameters
179
180 ! build overlap matrix
181 ALLOCATE (s(nterms, nterms))
182 s(:, :) = matmul(transpose(v_terms), v_terms)
183
184 ! diagonalize S
185 ALLOCATE (s_evals(nterms), s_evecs(nterms, nterms))
186 s_evecs(:, :) = s
187 CALL diamat_all(s_evecs, s_evals)
188
189 block_r = 0.0_dp
190 DO k = 1, nterms
191 v = pao%linpot_regu_delta/s_evals(k)
192 w = pao%linpot_regu_strength*min(1.0_dp, abs(v))
193 DO i = 1, nterms
194 DO j = 1, nterms
195 block_r(i, j) = block_r(i, j) + w*s_evecs(i, k)*s_evecs(j, k)
196 END DO
197 END DO
198 END DO
199
200 ! clean up
201 DEALLOCATE (s, s_evals, s_evecs)
202 END DO
203 CALL dbcsr_iterator_stop(iter)
204!$OMP END PARALLEL
205
206 CALL timestop(handle)
207 END SUBROUTINE pao_param_linpot_regularizer
208
209! **************************************************************************************************
210!> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
211!> \param pao ...
212! **************************************************************************************************
213 SUBROUTINE pao_param_linpot_preconditioner(pao)
214 TYPE(pao_env_type), POINTER :: pao
215
216 CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_linpot_preconditioner'
217
218 INTEGER :: acol, arow, handle, i, iatom, j, k, &
219 nterms
220 INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
221 LOGICAL :: found
222 REAL(dp) :: eval_capped
223 REAL(dp), ALLOCATABLE, DIMENSION(:) :: s_evals
224 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: s, s_evecs
225 REAL(dp), DIMENSION(:, :), POINTER :: block_precon, block_precon_inv, &
226 block_v_terms
227 TYPE(dbcsr_iterator_type) :: iter
228
229 CALL timeset(routinen, handle)
230
231 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot preconditioner"
232
233 CALL dbcsr_get_info(pao%matrix_V_terms, col_blk_size=blk_sizes_nterms)
234
235 CALL dbcsr_create(pao%matrix_precon, &
236 template=pao%matrix_V_terms, &
237 matrix_type="N", &
238 row_blk_size=blk_sizes_nterms, &
239 col_blk_size=blk_sizes_nterms, &
240 name="PAO matrix_precon")
241 CALL dbcsr_reserve_diag_blocks(pao%matrix_precon)
242
243 CALL dbcsr_create(pao%matrix_precon_inv, template=pao%matrix_precon, name="PAO matrix_precon_inv")
244 CALL dbcsr_reserve_diag_blocks(pao%matrix_precon_inv)
245
246!$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
247!$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)
248 CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
249 DO WHILE (dbcsr_iterator_blocks_left(iter))
250 CALL dbcsr_iterator_next_block(iter, arow, acol, block_v_terms)
251 iatom = arow; cpassert(arow == acol)
252 nterms = SIZE(block_v_terms, 2)
253 IF (nterms == 0) cycle ! protect against corner-case of zero pao parameters
254
255 CALL dbcsr_get_block_p(matrix=pao%matrix_precon, row=iatom, col=iatom, block=block_precon, found=found)
256 CALL dbcsr_get_block_p(matrix=pao%matrix_precon_inv, row=iatom, col=iatom, block=block_precon_inv, found=found)
257 cpassert(ASSOCIATED(block_precon))
258 cpassert(ASSOCIATED(block_precon_inv))
259
260 ALLOCATE (s(nterms, nterms))
261 s(:, :) = matmul(transpose(block_v_terms), block_v_terms)
262
263 ! diagonalize S
264 ALLOCATE (s_evals(nterms), s_evecs(nterms, nterms))
265 s_evecs(:, :) = s
266 CALL diamat_all(s_evecs, s_evals)
267
268 ! construct 1/Sqrt(S) and Sqrt(S)
269 block_precon = 0.0_dp
270 block_precon_inv = 0.0_dp
271 DO k = 1, nterms
272 eval_capped = max(pao%linpot_precon_delta, s_evals(k)) ! too small eigenvalues are hurtful
273 DO i = 1, nterms
274 DO j = 1, nterms
275 block_precon(i, j) = block_precon(i, j) + s_evecs(i, k)*s_evecs(j, k)/sqrt(eval_capped)
276 block_precon_inv(i, j) = block_precon_inv(i, j) + s_evecs(i, k)*s_evecs(j, k)*sqrt(eval_capped)
277 END DO
278 END DO
279 END DO
280
281 DEALLOCATE (s, s_evecs, s_evals)
282 END DO
283 CALL dbcsr_iterator_stop(iter)
284!$OMP END PARALLEL
285
286 CALL timestop(handle)
287 END SUBROUTINE pao_param_linpot_preconditioner
288
289! **************************************************************************************************
290!> \brief Finalize the linear potential parametrization
291!> \param pao ...
292! **************************************************************************************************
293 SUBROUTINE pao_param_finalize_linpot(pao)
294 TYPE(pao_env_type), POINTER :: pao
295
296 CALL dbcsr_release(pao%matrix_V_terms)
297 CALL dbcsr_release(pao%matrix_R)
298
299 IF (pao%precondition) THEN
300 CALL dbcsr_release(pao%matrix_precon)
301 CALL dbcsr_release(pao%matrix_precon_inv)
302 END IF
303
304 END SUBROUTINE pao_param_finalize_linpot
305
306! **************************************************************************************************
307!> \brief Returns the number of potential terms for given atomic kind
308!> \param pao ...
309!> \param qs_env ...
310!> \param ikind ...
311!> \param nparams ...
312! **************************************************************************************************
313 SUBROUTINE pao_param_count_linpot(pao, qs_env, ikind, nparams)
314 TYPE(pao_env_type), POINTER :: pao
315 TYPE(qs_environment_type), POINTER :: qs_env
316 INTEGER, INTENT(IN) :: ikind
317 INTEGER, INTENT(OUT) :: nparams
318
319 INTEGER :: pao_basis_size
320 TYPE(gto_basis_set_type), POINTER :: basis_set
321 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
322
323 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
324
325 CALL get_qs_kind(qs_kind_set(ikind), &
326 basis_set=basis_set, &
327 pao_basis_size=pao_basis_size)
328
329 IF (pao_basis_size == basis_set%nsgf) THEN
330 nparams = 0 ! pao disabled for iatom
331
332 ELSE
333 SELECT CASE (pao%parameterization)
334 CASE (pao_fock_param)
335 CALL linpot_full_count_terms(qs_env, ikind, nterms=nparams)
336 CASE (pao_rotinv_param)
337 CALL linpot_rotinv_count_terms(qs_env, ikind, nterms=nparams)
338 CASE DEFAULT
339 cpabort("unknown parameterization")
340 END SELECT
341 END IF
342
343 END SUBROUTINE pao_param_count_linpot
344
345! **************************************************************************************************
346!> \brief Takes current matrix_X and calculates the matrices A and B.
347!> \param pao ...
348!> \param qs_env ...
349!> \param ls_scf_env ...
350!> \param gradient ...
351!> \param penalty ...
352!> \param forces ...
353! **************************************************************************************************
354 SUBROUTINE pao_calc_ab_linpot(pao, qs_env, ls_scf_env, gradient, penalty, forces)
355 TYPE(pao_env_type), POINTER :: pao
356 TYPE(qs_environment_type), POINTER :: qs_env
357 TYPE(ls_scf_env_type), TARGET :: ls_scf_env
358 LOGICAL, INTENT(IN) :: gradient
359 REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
360 REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
361
362 CHARACTER(len=*), PARAMETER :: routinen = 'pao_calc_AB_linpot'
363
364 INTEGER :: handle
365 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
366 TYPE(dbcsr_type) :: matrix_m, matrix_u
367
368 CALL timeset(routinen, handle)
369 CALL get_qs_env(qs_env, matrix_s=matrix_s)
370 CALL dbcsr_create(matrix_u, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
371 CALL dbcsr_reserve_diag_blocks(matrix_u)
372
373 !TODO: move this condition into pao_calc_U, use matrix_N as template
374 IF (gradient) THEN
375 CALL pao_calc_grad_lnv_wrt_u(qs_env, ls_scf_env, matrix_m)
376 CALL pao_calc_u_linpot(pao, qs_env, matrix_u, matrix_m, pao%matrix_G, penalty, forces)
377 CALL dbcsr_release(matrix_m)
378 ELSE
379 CALL pao_calc_u_linpot(pao, qs_env, matrix_u, penalty=penalty)
380 END IF
381
382 CALL pao_calc_ab_from_u(pao, qs_env, ls_scf_env, matrix_u)
383 CALL dbcsr_release(matrix_u)
384 CALL timestop(handle)
385 END SUBROUTINE pao_calc_ab_linpot
386
387! **************************************************************************************************
388!> \brief Calculate new matrix U and optinally its gradient G
389!> \param pao ...
390!> \param qs_env ...
391!> \param matrix_U ...
392!> \param matrix_M ...
393!> \param matrix_G ...
394!> \param penalty ...
395!> \param forces ...
396! **************************************************************************************************
397 SUBROUTINE pao_calc_u_linpot(pao, qs_env, matrix_U, matrix_M, matrix_G, penalty, forces)
398 TYPE(pao_env_type), POINTER :: pao
399 TYPE(qs_environment_type), POINTER :: qs_env
400 TYPE(dbcsr_type) :: matrix_u
401 TYPE(dbcsr_type), OPTIONAL :: matrix_m, matrix_g
402 REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
403 REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
404
405 CHARACTER(len=*), PARAMETER :: routinen = 'pao_calc_U_linpot'
406
407 INTEGER :: acol, arow, group_handle, handle, iatom, &
408 kterm, n, natoms, nterms
409 LOGICAL :: found
410 REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
411 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: evals
412 REAL(dp), DIMENSION(:), POINTER :: vec_m2, vec_v
413 REAL(dp), DIMENSION(:, :), POINTER :: block_g, block_m1, block_m2, block_r, &
414 block_u, block_v, block_v_terms, &
415 block_x
416 REAL(dp), DIMENSION(:, :, :), POINTER :: m_blocks
417 REAL(kind=dp) :: regu_energy
418 TYPE(dbcsr_iterator_type) :: iter
419 TYPE(mp_comm_type) :: group
420
421 CALL timeset(routinen, handle)
422
423 cpassert(PRESENT(matrix_g) .EQV. PRESENT(matrix_m))
424
425 CALL get_qs_env(qs_env, natom=natoms)
426 ALLOCATE (gaps(natoms), evals(10, natoms)) ! printing 10 eigenvalues seems reasonable
427 evals(:, :) = 0.0_dp
428 gaps(:) = huge(1.0_dp)
429 regu_energy = 0.0_dp
430 CALL dbcsr_get_info(matrix_u, group=group_handle)
431 CALL group%set_handle(group_handle)
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
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
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:106
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:372
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:294
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:60
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:355
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:314
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:49
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:114
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
cp_control_types::dft_control_type
Definition cp_control_types.F:559
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:126
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171