(git:b279b6b)
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 ! **************************************************************************************************
12 MODULE pao_param_linpot
13  USE atomic_kind_types, ONLY: get_atomic_kind
14  USE basis_set_types, ONLY: gto_basis_set_type
15  USE cp_control_types, ONLY: dft_control_type
16  USE 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 
52 CONTAINS
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 
665 END 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