(git:e7e05ae)
pao_param_gth.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 Parametrization based on GTH pseudo potentials
10 !> \author Ole Schuett
11 ! **************************************************************************************************
12 MODULE pao_param_gth
13  USE arnoldi_api, ONLY: arnoldi_extremal
14  USE atomic_kind_types, ONLY: get_atomic_kind
15  USE basis_set_types, ONLY: gto_basis_set_type
16  USE cell_types, ONLY: cell_type,&
17  pbc
18  USE dbcsr_api, ONLY: &
19  dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
20  dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
21  dbcsr_p_type, dbcsr_release, dbcsr_reserve_all_blocks, dbcsr_reserve_diag_blocks, &
22  dbcsr_set, dbcsr_type
23  USE dm_ls_scf_types, ONLY: ls_scf_env_type
24  USE iterate_matrix, ONLY: matrix_sqrt_newton_schulz
25  USE kinds, ONLY: dp
26  USE machine, ONLY: m_flush
27  USE message_passing, ONLY: mp_comm_type
28  USE orbital_pointers, ONLY: init_orbital_pointers
29  USE pao_param_fock, ONLY: pao_calc_u_block_fock
30  USE pao_param_methods, ONLY: pao_calc_ab_from_u,&
31  pao_calc_grad_lnv_wrt_u
32  USE pao_potentials, ONLY: pao_calc_gaussian
33  USE pao_types, ONLY: pao_env_type
34  USE particle_types, ONLY: particle_type
35  USE qs_environment_types, ONLY: get_qs_env,&
36  qs_environment_type
37  USE qs_kind_types, ONLY: get_qs_kind,&
38  pao_potential_type,&
39  qs_kind_type
40 #include "./base/base_uses.f90"
41 
42  IMPLICIT NONE
43 
44  PRIVATE
45 
46  PUBLIC :: pao_param_init_gth, pao_param_finalize_gth, pao_calc_ab_gth
47  PUBLIC :: pao_param_count_gth, pao_param_initguess_gth
48 
49 CONTAINS
50 
51 ! **************************************************************************************************
52 !> \brief Initialize the linear potential parametrization
53 !> \param pao ...
54 !> \param qs_env ...
55 ! **************************************************************************************************
56  SUBROUTINE pao_param_init_gth(pao, qs_env)
57  TYPE(pao_env_type), POINTER :: pao
58  TYPE(qs_environment_type), POINTER :: qs_env
59 
60  CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_init_gth'
61 
62  INTEGER :: acol, arow, handle, iatom, idx, ikind, &
63  iterm, jatom, maxl, n, natoms
64  INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, nterms, &
65  row_blk_size
66  REAL(dp), DIMENSION(:, :), POINTER :: block_v_term, vec_v_terms
67  TYPE(dbcsr_iterator_type) :: iter
68  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
69  TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
70  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
71  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
72 
73  CALL timeset(routinen, handle)
74 
75  CALL get_qs_env(qs_env, &
76  natom=natoms, &
77  matrix_s=matrix_s, &
78  qs_kind_set=qs_kind_set, &
79  particle_set=particle_set)
80 
81  maxl = 0
82  ALLOCATE (row_blk_size(natoms), col_blk_size(natoms), nterms(natoms))
83  DO iatom = 1, natoms
84  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
85  CALL pao_param_count_gth(qs_env, ikind, nterms(iatom))
86  CALL get_qs_kind(qs_kind_set(ikind), pao_potentials=pao_potentials)
87  cpassert(SIZE(pao_potentials) == 1)
88  maxl = max(maxl, pao_potentials(1)%maxl)
89  END DO
90  CALL init_orbital_pointers(maxl) ! needs to be called before gth_calc_term()
91 
92  ! allocate matrix_V_terms
93  CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=blk_sizes_pri)
94  col_blk_size = sum(nterms)
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  CALL dbcsr_set(pao%matrix_V_terms, 0.0_dp)
104 
105  ! calculate and store poential terms
106 !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,blk_sizes_pri,natoms,nterms) &
107 !$OMP PRIVATE(iter,arow,acol,iatom,jatom,N,idx,vec_V_terms,block_V_term)
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, vec_v_terms)
111  iatom = arow; cpassert(arow == acol)
112  n = blk_sizes_pri(iatom)
113  DO jatom = 1, natoms
114  IF (jatom == iatom) cycle ! waste some storage to simplify things later
115  DO iterm = 1, nterms(jatom)
116  idx = sum(nterms(1:jatom - 1)) + iterm
117  block_v_term(1:n, 1:n) => vec_v_terms(:, idx) ! map column into matrix
118  CALL gth_calc_term(qs_env, block_v_term, iatom, jatom, iterm)
119  END DO
120  END DO
121  END DO
122  CALL dbcsr_iterator_stop(iter)
123 !$OMP END PARALLEL
124 
125  IF (pao%precondition) &
126  CALL pao_param_gth_preconditioner(pao, qs_env, nterms)
127 
128  DEALLOCATE (row_blk_size, col_blk_size, nterms)
129  CALL timestop(handle)
130  END SUBROUTINE pao_param_init_gth
131 
132 ! **************************************************************************************************
133 !> \brief Finalize the GTH potential parametrization
134 !> \param pao ...
135 ! **************************************************************************************************
136  SUBROUTINE pao_param_finalize_gth(pao)
137  TYPE(pao_env_type), POINTER :: pao
138 
139  CALL dbcsr_release(pao%matrix_V_terms)
140  IF (pao%precondition) THEN
141  CALL dbcsr_release(pao%matrix_precon)
142  CALL dbcsr_release(pao%matrix_precon_inv)
143  END IF
144 
145  END SUBROUTINE pao_param_finalize_gth
146 
147 ! **************************************************************************************************
148 !> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
149 !> \param pao ...
150 !> \param qs_env ...
151 !> \param nterms ...
152 ! **************************************************************************************************
153  SUBROUTINE pao_param_gth_preconditioner(pao, qs_env, nterms)
154  TYPE(pao_env_type), POINTER :: pao
155  TYPE(qs_environment_type), POINTER :: qs_env
156  INTEGER, DIMENSION(:), POINTER :: nterms
157 
158  CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_gth_preconditioner'
159 
160  INTEGER :: acol, arow, group_handle, handle, i, &
161  iatom, ioffset, j, jatom, joffset, m, &
162  n, natoms
163  LOGICAL :: arnoldi_converged, converged, found
164  REAL(dp) :: eval_max, eval_min
165  REAL(dp), DIMENSION(:, :), POINTER :: block, block_overlap, block_v_term
166  TYPE(dbcsr_iterator_type) :: iter
167  TYPE(dbcsr_type) :: matrix_gth_overlap
168  TYPE(ls_scf_env_type), POINTER :: ls_scf_env
169  TYPE(mp_comm_type) :: group
170 
171  CALL timeset(routinen, handle)
172 
173  CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env)
174  CALL dbcsr_get_info(pao%matrix_V_terms, group=group_handle)
175  CALL group%set_handle(group_handle)
176  natoms = SIZE(nterms)
177 
178  CALL dbcsr_create(matrix_gth_overlap, &
179  template=pao%matrix_V_terms, &
180  matrix_type="N", &
181  row_blk_size=nterms, &
182  col_blk_size=nterms)
183  CALL dbcsr_reserve_all_blocks(matrix_gth_overlap)
184  CALL dbcsr_set(matrix_gth_overlap, 0.0_dp)
185 
186  DO iatom = 1, natoms
187  DO jatom = 1, natoms
188  ioffset = sum(nterms(1:iatom - 1))
189  joffset = sum(nterms(1:jatom - 1))
190  n = nterms(iatom)
191  m = nterms(jatom)
192 
193  ALLOCATE (block(n, m))
194  block = 0.0_dp
195 
196  ! can't use OpenMP here block is a pointer and hence REDUCTION(+:block) does work
197  CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
198  DO WHILE (dbcsr_iterator_blocks_left(iter))
199  CALL dbcsr_iterator_next_block(iter, arow, acol, block_v_term)
200  cpassert(arow == acol)
201  DO i = 1, n
202  DO j = 1, m
203  block(i, j) = block(i, j) + sum(block_v_term(:, ioffset + i)*block_v_term(:, joffset + j))
204  END DO
205  END DO
206  END DO
207  CALL dbcsr_iterator_stop(iter)
208 
209  CALL group%sum(block)
210 
211  CALL dbcsr_get_block_p(matrix=matrix_gth_overlap, row=iatom, col=jatom, block=block_overlap, found=found)
212  IF (ASSOCIATED(block_overlap)) &
213  block_overlap = block
214 
215  DEALLOCATE (block)
216  END DO
217  END DO
218 
219  !TODO: good setting for arnoldi?
220  CALL arnoldi_extremal(matrix_gth_overlap, eval_max, eval_min, max_iter=100, &
221  threshold=1e-2_dp, converged=arnoldi_converged)
222  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| GTH-preconditioner converged, min, max, max/min:", &
223  arnoldi_converged, eval_min, eval_max, eval_max/eval_min
224 
225  CALL dbcsr_create(pao%matrix_precon, template=matrix_gth_overlap)
226  CALL dbcsr_create(pao%matrix_precon_inv, template=matrix_gth_overlap)
227 
228  CALL matrix_sqrt_newton_schulz(pao%matrix_precon_inv, pao%matrix_precon, matrix_gth_overlap, &
229  threshold=ls_scf_env%eps_filter, &
230  order=ls_scf_env%s_sqrt_order, &
231  max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
232  eps_lanczos=ls_scf_env%eps_lanczos, &
233  converged=converged)
234  CALL dbcsr_release(matrix_gth_overlap)
235 
236  IF (.NOT. converged) &
237  cpabort("PAO: Sqrt of GTH-preconditioner did not converge.")
238 
239  CALL timestop(handle)
240  END SUBROUTINE pao_param_gth_preconditioner
241 
242 ! **************************************************************************************************
243 !> \brief Takes current matrix_X and calculates the matrices A and B.
244 !> \param pao ...
245 !> \param qs_env ...
246 !> \param ls_scf_env ...
247 !> \param gradient ...
248 !> \param penalty ...
249 ! **************************************************************************************************
250  SUBROUTINE pao_calc_ab_gth(pao, qs_env, ls_scf_env, gradient, penalty)
251  TYPE(pao_env_type), POINTER :: pao
252  TYPE(qs_environment_type), POINTER :: qs_env
253  TYPE(ls_scf_env_type), TARGET :: ls_scf_env
254  LOGICAL, INTENT(IN) :: gradient
255  REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
256 
257  CHARACTER(len=*), PARAMETER :: routinen = 'pao_calc_AB_gth'
258 
259  INTEGER :: handle
260  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
261  TYPE(dbcsr_type) :: matrix_m, matrix_u
262 
263  CALL timeset(routinen, handle)
264  CALL get_qs_env(qs_env, matrix_s=matrix_s)
265  CALL dbcsr_create(matrix_u, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
266  CALL dbcsr_reserve_diag_blocks(matrix_u)
267 
268  !TODO: move this condition into pao_calc_U, use matrix_N as template
269  IF (gradient) THEN
270  CALL pao_calc_grad_lnv_wrt_u(qs_env, ls_scf_env, matrix_m)
271  CALL pao_calc_u_gth(pao, matrix_u, matrix_m, pao%matrix_G, penalty)
272  CALL dbcsr_release(matrix_m)
273  ELSE
274  CALL pao_calc_u_gth(pao, matrix_u, penalty=penalty)
275  END IF
276 
277  CALL pao_calc_ab_from_u(pao, qs_env, ls_scf_env, matrix_u)
278  CALL dbcsr_release(matrix_u)
279  CALL timestop(handle)
280  END SUBROUTINE pao_calc_ab_gth
281 
282 ! **************************************************************************************************
283 !> \brief Calculate new matrix U and optinally its gradient G
284 !> \param pao ...
285 !> \param matrix_U ...
286 !> \param matrix_M1 ...
287 !> \param matrix_G ...
288 !> \param penalty ...
289 ! **************************************************************************************************
290  SUBROUTINE pao_calc_u_gth(pao, matrix_U, matrix_M1, matrix_G, penalty)
291  TYPE(pao_env_type), POINTER :: pao
292  TYPE(dbcsr_type) :: matrix_u
293  TYPE(dbcsr_type), OPTIONAL :: matrix_m1, matrix_g
294  REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
295 
296  CHARACTER(len=*), PARAMETER :: routinen = 'pao_calc_U_gth'
297 
298  INTEGER :: acol, arow, group_handle, handle, iatom, &
299  idx, iterm, n, natoms
300  INTEGER, DIMENSION(:), POINTER :: nterms
301  LOGICAL :: found
302  REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
303  REAL(dp), DIMENSION(:), POINTER :: world_g, world_x
304  REAL(dp), DIMENSION(:, :), POINTER :: block_g, block_m1, block_m2, block_u, &
305  block_v, block_v_term, block_x, &
306  vec_v_terms
307  TYPE(dbcsr_iterator_type) :: iter
308  TYPE(mp_comm_type) :: group
309 
310  CALL timeset(routinen, handle)
311 
312  CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nterms, group=group_handle)
313  CALL group%set_handle(group_handle)
314  natoms = SIZE(nterms)
315  ALLOCATE (gaps(natoms))
316  gaps(:) = huge(dp)
317 
318  ! allocate arrays for world-view
319  ALLOCATE (world_x(sum(nterms)), world_g(sum(nterms)))
320  world_x = 0.0_dp; world_g = 0.0_dp
321 
322  ! collect world_X from atomic blocks
323  CALL dbcsr_iterator_start(iter, pao%matrix_X)
324  DO WHILE (dbcsr_iterator_blocks_left(iter))
325  CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
326  iatom = arow; cpassert(arow == acol)
327  idx = sum(nterms(1:iatom - 1))
328  world_x(idx + 1:idx + nterms(iatom)) = block_x(:, 1)
329  END DO
330  CALL dbcsr_iterator_stop(iter)
331  CALL group%sum(world_x) ! sync world view across MPI ranks
332 
333  ! loop over atoms
334  CALL dbcsr_iterator_start(iter, matrix_u)
335  DO WHILE (dbcsr_iterator_blocks_left(iter))
336  CALL dbcsr_iterator_next_block(iter, arow, acol, block_u)
337  iatom = arow; cpassert(arow == acol)
338  n = SIZE(block_u, 1)
339  CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=vec_v_terms, found=found)
340  cpassert(ASSOCIATED(vec_v_terms))
341 
342  ! calculate potential V of i'th atom
343  ALLOCATE (block_v(n, n))
344  block_v = 0.0_dp
345  DO iterm = 1, SIZE(world_x)
346  block_v_term(1:n, 1:n) => vec_v_terms(:, iterm) ! map column into matrix
347  block_v = block_v + world_x(iterm)*block_v_term
348  END DO
349 
350  ! calculate gradient block of i'th atom
351  IF (.NOT. PRESENT(matrix_g)) THEN
352  CALL pao_calc_u_block_fock(pao, iatom=iatom, penalty=penalty, v=block_v, u=block_u, gap=gaps(iatom))
353 
354  ELSE ! TURNING POINT (if calc grad) ------------------------------------
355  cpassert(PRESENT(matrix_m1))
356  CALL dbcsr_get_block_p(matrix=matrix_m1, row=iatom, col=iatom, block=block_m1, found=found)
357  ALLOCATE (block_m2(n, n))
358  CALL pao_calc_u_block_fock(pao, iatom=iatom, penalty=penalty, v=block_v, u=block_u, &
359  m1=block_m1, g=block_m2, gap=gaps(iatom))
360  DO iterm = 1, SIZE(world_g)
361  block_v_term(1:n, 1:n) => vec_v_terms(:, iterm) ! map column into matrix
362  world_g(iterm) = world_g(iterm) + sum(block_v_term*block_m2)
363  END DO
364  DEALLOCATE (block_m2)
365  END IF
366  DEALLOCATE (block_v)
367  END DO
368  CALL dbcsr_iterator_stop(iter)
369 
370  ! distribute world_G across atomic blocks
371  IF (PRESENT(matrix_g)) THEN
372  CALL group%sum(world_g) ! sync world view across MPI ranks
373  CALL dbcsr_iterator_start(iter, matrix_g)
374  DO WHILE (dbcsr_iterator_blocks_left(iter))
375  CALL dbcsr_iterator_next_block(iter, arow, acol, block_g)
376  iatom = arow; cpassert(arow == acol)
377  idx = sum(nterms(1:iatom - 1))
378  block_g(:, 1) = world_g(idx + 1:idx + nterms(iatom))
379  END DO
380  CALL dbcsr_iterator_stop(iter)
381  END IF
382 
383  DEALLOCATE (world_x, world_g)
384 
385  ! sum penalty energies across ranks
386  IF (PRESENT(penalty)) &
387  CALL group%sum(penalty)
388 
389  ! print homo-lumo gap encountered by fock-layer
390  CALL group%min(gaps)
391  IF (pao%iw_gap > 0) THEN
392  DO iatom = 1, natoms
393  WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
394  END DO
395  CALL m_flush(pao%iw_gap)
396  END IF
397 
398  ! one-line summary
399  IF (pao%iw > 0) THEN
400  WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", minval(gaps), " for atom:", minloc(gaps)
401  END IF
402 
403  DEALLOCATE (gaps)
404  CALL timestop(handle)
405 
406  END SUBROUTINE pao_calc_u_gth
407 
408 ! **************************************************************************************************
409 !> \brief Returns the number of parameters for given atomic kind
410 !> \param qs_env ...
411 !> \param ikind ...
412 !> \param nparams ...
413 ! **************************************************************************************************
414  SUBROUTINE pao_param_count_gth(qs_env, ikind, nparams)
415  TYPE(qs_environment_type), POINTER :: qs_env
416  INTEGER, INTENT(IN) :: ikind
417  INTEGER, INTENT(OUT) :: nparams
418 
419  INTEGER :: max_projector, maxl, ncombis
420  TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
421  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
422 
423  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
424  CALL get_qs_kind(qs_kind_set(ikind), pao_potentials=pao_potentials)
425 
426  IF (SIZE(pao_potentials) /= 1) &
427  cpabort("GTH parametrization requires exactly one PAO_POTENTIAL section per KIND")
428 
429  max_projector = pao_potentials(1)%max_projector
430  maxl = pao_potentials(1)%maxl
431 
432  IF (maxl < 0) &
433  cpabort("GTH parametrization requires non-negative PAO_POTENTIAL%MAXL")
434 
435  IF (max_projector < 0) &
436  cpabort("GTH parametrization requires non-negative PAO_POTENTIAL%MAX_PROJECTOR")
437 
438  IF (mod(maxl, 2) /= 0) &
439  cpabort("GTH parametrization requires even-numbered PAO_POTENTIAL%MAXL")
440 
441  ncombis = (max_projector + 1)*(max_projector + 2)/2
442  nparams = ncombis*(maxl/2 + 1)
443 
444  END SUBROUTINE pao_param_count_gth
445 
446 ! **************************************************************************************************
447 !> \brief Fills the given block_V with the requested potential term
448 !> \param qs_env ...
449 !> \param block_V ...
450 !> \param iatom ...
451 !> \param jatom ...
452 !> \param kterm ...
453 ! **************************************************************************************************
454  SUBROUTINE gth_calc_term(qs_env, block_V, iatom, jatom, kterm)
455  TYPE(qs_environment_type), POINTER :: qs_env
456  REAL(dp), DIMENSION(:, :), INTENT(OUT) :: block_v
457  INTEGER, INTENT(IN) :: iatom, jatom, kterm
458 
459  INTEGER :: c, ikind, jkind, lpot, max_l, min_l, &
460  pot_max_projector, pot_maxl
461  REAL(dp), DIMENSION(3) :: ra, rab, rb
462  REAL(kind=dp) :: pot_beta
463  TYPE(cell_type), POINTER :: cell
464  TYPE(gto_basis_set_type), POINTER :: basis_set
465  TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
466  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
467  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
468 
469  CALL get_qs_env(qs_env, &
470  cell=cell, &
471  particle_set=particle_set, &
472  qs_kind_set=qs_kind_set)
473 
474  ! get GTH-settings from remote atom
475  CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
476  CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=pao_potentials)
477  cpassert(SIZE(pao_potentials) == 1)
478  pot_max_projector = pao_potentials(1)%max_projector
479  pot_maxl = pao_potentials(1)%maxl
480  pot_beta = pao_potentials(1)%beta
481 
482  c = 0
483  outer: DO lpot = 0, pot_maxl, 2
484  DO max_l = 0, pot_max_projector
485  DO min_l = 0, max_l
486  c = c + 1
487  IF (c == kterm) EXIT outer
488  END DO
489  END DO
490  END DO outer
491 
492  ! get basis-set of central atom
493  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
494  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set)
495 
496  ra = particle_set(iatom)%r
497  rb = particle_set(jatom)%r
498  rab = pbc(ra, rb, cell)
499 
500  block_v = 0.0_dp
501  CALL pao_calc_gaussian(basis_set, block_v, rab=rab, lpot=lpot, &
502  min_l=min_l, max_l=max_l, beta=pot_beta, weight=1.0_dp)
503 
504  END SUBROUTINE gth_calc_term
505 
506 ! **************************************************************************************************
507 !> \brief Calculate initial guess for matrix_X
508 !> \param pao ...
509 ! **************************************************************************************************
510  SUBROUTINE pao_param_initguess_gth(pao)
511  TYPE(pao_env_type), POINTER :: pao
512 
513  INTEGER :: acol, arow
514  REAL(dp), DIMENSION(:, :), POINTER :: block_x
515  TYPE(dbcsr_iterator_type) :: iter
516 
517 !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
518 !$OMP PRIVATE(iter,arow,acol,block_X)
519  CALL dbcsr_iterator_start(iter, pao%matrix_X)
520  DO WHILE (dbcsr_iterator_blocks_left(iter))
521  CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
522  cpassert(arow == acol)
523  cpassert(SIZE(block_x, 2) == 1)
524 
525  ! a simplistic guess, which at least makes the atom visible to others
526  block_x = 0.0_dp
527  block_x(1, 1) = 0.01_dp
528  END DO
529  CALL dbcsr_iterator_stop(iter)
530 !$OMP END PARALLEL
531 
532  END SUBROUTINE pao_param_initguess_gth
533 
534 END MODULE pao_param_gth
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition: grid_common.h:153
arnoldi_api
arnoldi iteration using dbcsr
Definition: arnoldi_api.F:15
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition: arnoldi_api.F:254
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
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
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
iterate_matrix
Routines useful for iterative matrix calculations.
Definition: iterate_matrix.F:13
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition: iterate_matrix.F:1624
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
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
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_gth
Parametrization based on GTH pseudo potentials.
Definition: pao_param_gth.F:12
pao_param_gth::pao_calc_ab_gth
subroutine, public pao_calc_ab_gth(pao, qs_env, ls_scf_env, gradient, penalty)
Takes current matrix_X and calculates the matrices A and B.
Definition: pao_param_gth.F:251
pao_param_gth::pao_param_initguess_gth
subroutine, public pao_param_initguess_gth(pao)
Calculate initial guess for matrix_X.
Definition: pao_param_gth.F:511
pao_param_gth::pao_param_count_gth
subroutine, public pao_param_count_gth(qs_env, ikind, nparams)
Returns the number of parameters for given atomic kind.
Definition: pao_param_gth.F:415
pao_param_gth::pao_param_init_gth
subroutine, public pao_param_init_gth(pao, qs_env)
Initialize the linear potential parametrization.
Definition: pao_param_gth.F:57
pao_param_gth::pao_param_finalize_gth
subroutine, public pao_param_finalize_gth(pao)
Finalize the GTH potential parametrization.
Definition: pao_param_gth.F:137
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_calc_gaussian
subroutine, public pao_calc_gaussian(basis_set, block_V, block_D, Rab, lpot, beta, weight, min_shell, max_shell, min_l, max_l)
Calculates potential term of the form r**lpot * Exp(-beta*r**2) One needs to call init_orbital_pointe...
Definition: pao_potentials.F:102
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