Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 cp_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_set, dbcsr_type
22 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks,&
23 : dbcsr_reserve_diag_blocks
24 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
25 : USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
26 : USE kinds, ONLY: dp
27 : USE machine, ONLY: m_flush
28 : USE message_passing, ONLY: mp_comm_type
29 : USE orbital_pointers, ONLY: init_orbital_pointers
30 : USE pao_param_fock, ONLY: pao_calc_U_block_fock
31 : USE pao_param_methods, ONLY: pao_calc_AB_from_U,&
32 : pao_calc_grad_lnv_wrt_U
33 : USE pao_potentials, ONLY: pao_calc_gaussian
34 : USE pao_types, ONLY: pao_env_type
35 : USE particle_types, ONLY: particle_type
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type
38 : USE qs_kind_types, ONLY: get_qs_kind,&
39 : pao_potential_type,&
40 : qs_kind_type
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : PUBLIC :: pao_param_init_gth, pao_param_finalize_gth, pao_calc_AB_gth
48 : PUBLIC :: pao_param_count_gth, pao_param_initguess_gth
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Initialize the linear potential parametrization
54 : !> \param pao ...
55 : !> \param qs_env ...
56 : ! **************************************************************************************************
57 10 : SUBROUTINE pao_param_init_gth(pao, qs_env)
58 : TYPE(pao_env_type), POINTER :: pao
59 : TYPE(qs_environment_type), POINTER :: qs_env
60 :
61 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_init_gth'
62 :
63 : INTEGER :: acol, arow, handle, iatom, idx, ikind, &
64 : iterm, jatom, maxl, n, natoms
65 10 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, nterms, &
66 10 : row_blk_size
67 10 : REAL(dp), DIMENSION(:, :), POINTER :: block_V_term, vec_V_terms
68 : TYPE(dbcsr_iterator_type) :: iter
69 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
70 10 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
71 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
72 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
73 :
74 10 : CALL timeset(routineN, handle)
75 :
76 : CALL get_qs_env(qs_env, &
77 : natom=natoms, &
78 : matrix_s=matrix_s, &
79 : qs_kind_set=qs_kind_set, &
80 10 : particle_set=particle_set)
81 :
82 10 : maxl = 0
83 50 : ALLOCATE (row_blk_size(natoms), col_blk_size(natoms), nterms(natoms))
84 32 : DO iatom = 1, natoms
85 22 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
86 22 : CALL pao_param_count_gth(qs_env, ikind, nterms(iatom))
87 22 : CALL get_qs_kind(qs_kind_set(ikind), pao_potentials=pao_potentials)
88 22 : CPASSERT(SIZE(pao_potentials) == 1)
89 54 : maxl = MAX(maxl, pao_potentials(1)%maxl)
90 : END DO
91 10 : CALL init_orbital_pointers(maxl) ! needs to be called before gth_calc_term()
92 :
93 : ! allocate matrix_V_terms
94 10 : CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=blk_sizes_pri)
95 54 : col_blk_size = SUM(nterms)
96 64 : row_blk_size = blk_sizes_pri**2
97 : CALL dbcsr_create(pao%matrix_V_terms, &
98 : name="PAO matrix_V_terms", &
99 : dist=pao%diag_distribution, &
100 : matrix_type="N", &
101 : row_blk_size=row_blk_size, &
102 10 : col_blk_size=col_blk_size)
103 10 : CALL dbcsr_reserve_diag_blocks(pao%matrix_V_terms)
104 10 : CALL dbcsr_set(pao%matrix_V_terms, 0.0_dp)
105 :
106 : ! calculate and store poential terms
107 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,blk_sizes_pri,natoms,nterms) &
108 10 : !$OMP PRIVATE(iter,arow,acol,iatom,jatom,N,idx,vec_V_terms,block_V_term)
109 : CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
110 : DO WHILE (dbcsr_iterator_blocks_left(iter))
111 : CALL dbcsr_iterator_next_block(iter, arow, acol, vec_V_terms)
112 : iatom = arow; CPASSERT(arow == acol)
113 : n = blk_sizes_pri(iatom)
114 : DO jatom = 1, natoms
115 : IF (jatom == iatom) CYCLE ! waste some storage to simplify things later
116 : DO iterm = 1, nterms(jatom)
117 : idx = SUM(nterms(1:jatom - 1)) + iterm
118 : block_V_term(1:n, 1:n) => vec_V_terms(:, idx) ! map column into matrix
119 : CALL gth_calc_term(qs_env, block_V_term, iatom, jatom, iterm)
120 : END DO
121 : END DO
122 : END DO
123 : CALL dbcsr_iterator_stop(iter)
124 : !$OMP END PARALLEL
125 :
126 10 : IF (pao%precondition) THEN
127 4 : CALL pao_param_gth_preconditioner(pao, qs_env, nterms)
128 : END IF
129 :
130 10 : DEALLOCATE (row_blk_size, col_blk_size, nterms)
131 10 : CALL timestop(handle)
132 10 : END SUBROUTINE pao_param_init_gth
133 :
134 : ! **************************************************************************************************
135 : !> \brief Finalize the GTH potential parametrization
136 : !> \param pao ...
137 : ! **************************************************************************************************
138 10 : SUBROUTINE pao_param_finalize_gth(pao)
139 : TYPE(pao_env_type), POINTER :: pao
140 :
141 10 : CALL dbcsr_release(pao%matrix_V_terms)
142 10 : IF (pao%precondition) THEN
143 4 : CALL dbcsr_release(pao%matrix_precon)
144 4 : CALL dbcsr_release(pao%matrix_precon_inv)
145 : END IF
146 :
147 10 : END SUBROUTINE pao_param_finalize_gth
148 :
149 : ! **************************************************************************************************
150 : !> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
151 : !> \param pao ...
152 : !> \param qs_env ...
153 : !> \param nterms ...
154 : ! **************************************************************************************************
155 8 : SUBROUTINE pao_param_gth_preconditioner(pao, qs_env, nterms)
156 : TYPE(pao_env_type), POINTER :: pao
157 : TYPE(qs_environment_type), POINTER :: qs_env
158 : INTEGER, DIMENSION(:), POINTER :: nterms
159 :
160 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_gth_preconditioner'
161 :
162 : INTEGER :: acol, arow, handle, i, iatom, ioffset, &
163 : j, jatom, joffset, m, n, natoms
164 : LOGICAL :: arnoldi_converged, converged, found
165 : REAL(dp) :: eval_max, eval_min
166 4 : REAL(dp), DIMENSION(:, :), POINTER :: block, block_overlap, block_V_term
167 : TYPE(dbcsr_iterator_type) :: iter
168 : TYPE(dbcsr_type) :: matrix_gth_overlap
169 : TYPE(ls_scf_env_type), POINTER :: ls_scf_env
170 : TYPE(mp_comm_type) :: group
171 :
172 4 : CALL timeset(routineN, handle)
173 :
174 4 : CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env)
175 4 : CALL dbcsr_get_info(pao%matrix_V_terms, group=group)
176 4 : 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 4 : col_blk_size=nterms)
183 4 : CALL dbcsr_reserve_all_blocks(matrix_gth_overlap)
184 4 : CALL dbcsr_set(matrix_gth_overlap, 0.0_dp)
185 :
186 16 : DO iatom = 1, natoms
187 52 : DO jatom = 1, natoms
188 72 : ioffset = SUM(nterms(1:iatom - 1))
189 72 : joffset = SUM(nterms(1:jatom - 1))
190 36 : n = nterms(iatom)
191 36 : m = nterms(jatom)
192 :
193 144 : ALLOCATE (block(n, m))
194 3996 : block = 0.0_dp
195 :
196 : ! can't use OpenMP here block is a pointer and hence REDUCTION(+:block) does work
197 36 : CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
198 90 : DO WHILE (dbcsr_iterator_blocks_left(iter))
199 54 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_V_term)
200 54 : CPASSERT(arow == acol)
201 630 : DO i = 1, n
202 5994 : DO j = 1, m
203 400140 : 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 36 : CALL dbcsr_iterator_stop(iter)
208 :
209 7956 : CALL group%sum(block)
210 :
211 36 : CALL dbcsr_get_block_p(matrix=matrix_gth_overlap, row=iatom, col=jatom, block=block_overlap, found=found)
212 36 : IF (ASSOCIATED(block_overlap)) THEN
213 3996 : block_overlap = block
214 : END IF
215 :
216 120 : DEALLOCATE (block)
217 : END DO
218 : END DO
219 :
220 : !TODO: good setting for arnoldi?
221 : CALL arnoldi_extremal(matrix_gth_overlap, eval_max, eval_min, max_iter=100, &
222 4 : threshold=1e-2_dp, converged=arnoldi_converged)
223 6 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| GTH-preconditioner converged, min, max, max/min:", &
224 4 : arnoldi_converged, eval_min, eval_max, eval_max/eval_min
225 :
226 4 : CALL dbcsr_create(pao%matrix_precon, template=matrix_gth_overlap)
227 4 : CALL dbcsr_create(pao%matrix_precon_inv, template=matrix_gth_overlap)
228 :
229 : CALL matrix_sqrt_Newton_Schulz(pao%matrix_precon_inv, pao%matrix_precon, matrix_gth_overlap, &
230 : threshold=ls_scf_env%eps_filter, &
231 : order=ls_scf_env%s_sqrt_order, &
232 : max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
233 : eps_lanczos=ls_scf_env%eps_lanczos, &
234 4 : converged=converged)
235 4 : CALL dbcsr_release(matrix_gth_overlap)
236 :
237 4 : IF (.NOT. converged) THEN
238 0 : CPABORT("PAO: Sqrt of GTH-preconditioner did not converge.")
239 : END IF
240 :
241 4 : CALL timestop(handle)
242 4 : END SUBROUTINE pao_param_gth_preconditioner
243 :
244 : ! **************************************************************************************************
245 : !> \brief Takes current matrix_X and calculates the matrices A and B.
246 : !> \param pao ...
247 : !> \param qs_env ...
248 : !> \param ls_scf_env ...
249 : !> \param gradient ...
250 : !> \param penalty ...
251 : ! **************************************************************************************************
252 2152 : SUBROUTINE pao_calc_AB_gth(pao, qs_env, ls_scf_env, gradient, penalty)
253 : TYPE(pao_env_type), POINTER :: pao
254 : TYPE(qs_environment_type), POINTER :: qs_env
255 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
256 : LOGICAL, INTENT(IN) :: gradient
257 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
258 :
259 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_AB_gth'
260 :
261 : INTEGER :: handle
262 2152 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
263 : TYPE(dbcsr_type) :: matrix_M, matrix_U
264 :
265 2152 : CALL timeset(routineN, handle)
266 2152 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
267 2152 : CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
268 2152 : CALL dbcsr_reserve_diag_blocks(matrix_U)
269 :
270 : !TODO: move this condition into pao_calc_U, use matrix_N as template
271 2152 : IF (gradient) THEN
272 322 : CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
273 322 : CALL pao_calc_U_gth(pao, matrix_U, matrix_M, pao%matrix_G, penalty)
274 322 : CALL dbcsr_release(matrix_M)
275 : ELSE
276 1830 : CALL pao_calc_U_gth(pao, matrix_U, penalty=penalty)
277 : END IF
278 :
279 2152 : CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
280 2152 : CALL dbcsr_release(matrix_U)
281 2152 : CALL timestop(handle)
282 2152 : END SUBROUTINE pao_calc_AB_gth
283 :
284 : ! **************************************************************************************************
285 : !> \brief Calculate new matrix U and optinally its gradient G
286 : !> \param pao ...
287 : !> \param matrix_U ...
288 : !> \param matrix_M1 ...
289 : !> \param matrix_G ...
290 : !> \param penalty ...
291 : ! **************************************************************************************************
292 2152 : SUBROUTINE pao_calc_U_gth(pao, matrix_U, matrix_M1, matrix_G, penalty)
293 : TYPE(pao_env_type), POINTER :: pao
294 : TYPE(dbcsr_type) :: matrix_U
295 : TYPE(dbcsr_type), OPTIONAL :: matrix_M1, matrix_G
296 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
297 :
298 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_U_gth'
299 :
300 : INTEGER :: acol, arow, handle, iatom, idx, iterm, &
301 : n, natoms
302 2152 : INTEGER, DIMENSION(:), POINTER :: nterms
303 : LOGICAL :: found
304 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
305 2152 : REAL(dp), DIMENSION(:), POINTER :: world_G, world_X
306 2152 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_M1, block_M2, block_U, &
307 2152 : block_V, block_V_term, block_X, &
308 2152 : vec_V_terms
309 : TYPE(dbcsr_iterator_type) :: iter
310 : TYPE(mp_comm_type) :: group
311 :
312 2152 : CALL timeset(routineN, handle)
313 :
314 2152 : CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nterms, group=group)
315 2152 : natoms = SIZE(nterms)
316 6456 : ALLOCATE (gaps(natoms))
317 7620 : gaps(:) = HUGE(dp)
318 :
319 : ! allocate arrays for world-view
320 23848 : ALLOCATE (world_X(SUM(nterms)), world_G(SUM(nterms)))
321 204320 : world_X = 0.0_dp; world_G = 0.0_dp
322 :
323 : ! collect world_X from atomic blocks
324 2152 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
325 4886 : DO WHILE (dbcsr_iterator_blocks_left(iter))
326 2734 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
327 2734 : iatom = arow; CPASSERT(arow == acol)
328 4975 : idx = SUM(nterms(1:iatom - 1))
329 107628 : world_X(idx + 1:idx + nterms(iatom)) = block_X(:, 1)
330 : END DO
331 2152 : CALL dbcsr_iterator_stop(iter)
332 202168 : CALL group%sum(world_X) ! sync world view across MPI ranks
333 :
334 : ! loop over atoms
335 2152 : CALL dbcsr_iterator_start(iter, matrix_U)
336 4886 : DO WHILE (dbcsr_iterator_blocks_left(iter))
337 2734 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_U)
338 2734 : iatom = arow; CPASSERT(arow == acol)
339 2734 : n = SIZE(block_U, 1)
340 2734 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=vec_V_terms, found=found)
341 2734 : CPASSERT(ASSOCIATED(vec_V_terms))
342 :
343 : ! calculate potential V of i'th atom
344 10936 : ALLOCATE (block_V(n, n))
345 173370 : block_V = 0.0_dp
346 120226 : DO iterm = 1, SIZE(world_X)
347 117492 : block_V_term(1:n, 1:n) => vec_V_terms(:, iterm) ! map column into matrix
348 12604198 : block_V = block_V + world_X(iterm)*block_V_term
349 : END DO
350 :
351 : ! calculate gradient block of i'th atom
352 2734 : IF (.NOT. PRESENT(matrix_G)) THEN
353 2288 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, gap=gaps(iatom))
354 :
355 : ELSE ! TURNING POINT (if calc grad) ------------------------------------
356 446 : CPASSERT(PRESENT(matrix_M1))
357 446 : CALL dbcsr_get_block_p(matrix=matrix_M1, row=iatom, col=iatom, block=block_M1, found=found)
358 1338 : ALLOCATE (block_M2(n, n))
359 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
360 446 : M1=block_M1, G=block_M2, gap=gaps(iatom))
361 16910 : DO iterm = 1, SIZE(world_G)
362 16464 : block_V_term(1:n, 1:n) => vec_V_terms(:, iterm) ! map column into matrix
363 1076270 : world_G(iterm) = world_G(iterm) + SUM(block_V_term*block_M2)
364 : END DO
365 892 : DEALLOCATE (block_M2)
366 : END IF
367 10354 : DEALLOCATE (block_V)
368 : END DO
369 2152 : CALL dbcsr_iterator_stop(iter)
370 :
371 : ! distribute world_G across atomic blocks
372 2152 : IF (PRESENT(matrix_G)) THEN
373 25810 : CALL group%sum(world_G) ! sync world view across MPI ranks
374 322 : CALL dbcsr_iterator_start(iter, matrix_G)
375 768 : DO WHILE (dbcsr_iterator_blocks_left(iter))
376 446 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_G)
377 446 : iatom = arow; CPASSERT(arow == acol)
378 855 : idx = SUM(nterms(1:iatom - 1))
379 13958 : block_G(:, 1) = world_G(idx + 1:idx + nterms(iatom))
380 : END DO
381 322 : CALL dbcsr_iterator_stop(iter)
382 : END IF
383 :
384 2152 : DEALLOCATE (world_X, world_G)
385 :
386 : ! sum penalty energies across ranks
387 2152 : IF (PRESENT(penalty)) THEN
388 2142 : CALL group%sum(penalty)
389 : END IF
390 :
391 : ! print homo-lumo gap encountered by fock-layer
392 2152 : CALL group%min(gaps)
393 2152 : IF (pao%iw_gap > 0) THEN
394 2208 : DO iatom = 1, natoms
395 2208 : WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
396 : END DO
397 552 : CALL m_flush(pao%iw_gap)
398 : END IF
399 :
400 : ! one-line summary
401 2152 : IF (pao%iw > 0) THEN
402 7620 : WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", MINVAL(gaps), " for atom:", MINLOC(gaps)
403 : END IF
404 :
405 2152 : DEALLOCATE (gaps)
406 2152 : CALL timestop(handle)
407 :
408 6456 : END SUBROUTINE pao_calc_U_gth
409 :
410 : ! **************************************************************************************************
411 : !> \brief Returns the number of parameters for given atomic kind
412 : !> \param qs_env ...
413 : !> \param ikind ...
414 : !> \param nparams ...
415 : ! **************************************************************************************************
416 44 : SUBROUTINE pao_param_count_gth(qs_env, ikind, nparams)
417 : TYPE(qs_environment_type), POINTER :: qs_env
418 : INTEGER, INTENT(IN) :: ikind
419 : INTEGER, INTENT(OUT) :: nparams
420 :
421 : INTEGER :: max_projector, maxl, ncombis
422 44 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
423 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
424 :
425 44 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
426 44 : CALL get_qs_kind(qs_kind_set(ikind), pao_potentials=pao_potentials)
427 :
428 44 : IF (SIZE(pao_potentials) /= 1) THEN
429 0 : CPABORT("GTH parametrization requires exactly one PAO_POTENTIAL section per KIND")
430 : END IF
431 :
432 44 : max_projector = pao_potentials(1)%max_projector
433 44 : maxl = pao_potentials(1)%maxl
434 :
435 44 : IF (maxl < 0) THEN
436 0 : CPABORT("GTH parametrization requires non-negative PAO_POTENTIAL%MAXL")
437 : END IF
438 :
439 44 : IF (max_projector < 0) THEN
440 0 : CPABORT("GTH parametrization requires non-negative PAO_POTENTIAL%MAX_PROJECTOR")
441 : END IF
442 :
443 44 : IF (MOD(maxl, 2) /= 0) THEN
444 0 : CPABORT("GTH parametrization requires even-numbered PAO_POTENTIAL%MAXL")
445 : END IF
446 :
447 44 : ncombis = (max_projector + 1)*(max_projector + 2)/2
448 44 : nparams = ncombis*(maxl/2 + 1)
449 :
450 44 : END SUBROUTINE pao_param_count_gth
451 :
452 : ! **************************************************************************************************
453 : !> \brief Fills the given block_V with the requested potential term
454 : !> \param qs_env ...
455 : !> \param block_V ...
456 : !> \param iatom ...
457 : !> \param jatom ...
458 : !> \param kterm ...
459 : ! **************************************************************************************************
460 252 : SUBROUTINE gth_calc_term(qs_env, block_V, iatom, jatom, kterm)
461 : TYPE(qs_environment_type), POINTER :: qs_env
462 : REAL(dp), DIMENSION(:, :), INTENT(OUT) :: block_V
463 : INTEGER, INTENT(IN) :: iatom, jatom, kterm
464 :
465 : INTEGER :: c, ikind, jkind, lpot, max_l, min_l, &
466 : pot_max_projector, pot_maxl
467 : REAL(dp), DIMENSION(3) :: Ra, Rab, Rb
468 : REAL(KIND=dp) :: pot_beta
469 : TYPE(cell_type), POINTER :: cell
470 : TYPE(gto_basis_set_type), POINTER :: basis_set
471 252 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
472 252 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
473 252 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
474 :
475 : CALL get_qs_env(qs_env, &
476 : cell=cell, &
477 : particle_set=particle_set, &
478 252 : qs_kind_set=qs_kind_set)
479 :
480 : ! get GTH-settings from remote atom
481 252 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
482 252 : CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=pao_potentials)
483 252 : CPASSERT(SIZE(pao_potentials) == 1)
484 252 : pot_max_projector = pao_potentials(1)%max_projector
485 252 : pot_maxl = pao_potentials(1)%maxl
486 252 : pot_beta = pao_potentials(1)%beta
487 :
488 252 : c = 0
489 612 : outer: DO lpot = 0, pot_maxl, 2
490 2252 : DO max_l = 0, pot_max_projector
491 4718 : DO min_l = 0, max_l
492 2970 : c = c + 1
493 4358 : IF (c == kterm) EXIT outer
494 : END DO
495 : END DO
496 : END DO outer
497 :
498 : ! get basis-set of central atom
499 252 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
500 252 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set)
501 :
502 1008 : Ra = particle_set(iatom)%r
503 1008 : Rb = particle_set(jatom)%r
504 252 : Rab = pbc(ra, rb, cell)
505 :
506 15108 : block_V = 0.0_dp
507 : CALL pao_calc_gaussian(basis_set, block_V, Rab=Rab, lpot=lpot, &
508 252 : min_l=min_l, max_l=max_l, beta=pot_beta, weight=1.0_dp)
509 :
510 252 : END SUBROUTINE gth_calc_term
511 :
512 : ! **************************************************************************************************
513 : !> \brief Calculate initial guess for matrix_X
514 : !> \param pao ...
515 : ! **************************************************************************************************
516 10 : SUBROUTINE pao_param_initguess_gth(pao)
517 : TYPE(pao_env_type), POINTER :: pao
518 :
519 : INTEGER :: acol, arow
520 10 : REAL(dp), DIMENSION(:, :), POINTER :: block_X
521 : TYPE(dbcsr_iterator_type) :: iter
522 :
523 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
524 10 : !$OMP PRIVATE(iter,arow,acol,block_X)
525 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
526 : DO WHILE (dbcsr_iterator_blocks_left(iter))
527 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
528 : CPASSERT(arow == acol)
529 : CPASSERT(SIZE(block_X, 2) == 1)
530 :
531 : ! a simplistic guess, which at least makes the atom visible to others
532 : block_X = 0.0_dp
533 : block_X(1, 1) = 0.01_dp
534 : END DO
535 : CALL dbcsr_iterator_stop(iter)
536 : !$OMP END PARALLEL
537 :
538 10 : END SUBROUTINE pao_param_initguess_gth
539 :
540 : END MODULE pao_param_gth
|
