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