Line data Source code
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 10 : 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 10 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, nterms, &
65 10 : row_blk_size
66 10 : REAL(dp), DIMENSION(:, :), POINTER :: block_V_term, vec_V_terms
67 : TYPE(dbcsr_iterator_type) :: iter
68 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
69 10 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
70 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
71 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
72 :
73 10 : 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 10 : particle_set=particle_set)
80 :
81 10 : maxl = 0
82 50 : ALLOCATE (row_blk_size(natoms), col_blk_size(natoms), nterms(natoms))
83 32 : DO iatom = 1, natoms
84 22 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
85 22 : CALL pao_param_count_gth(qs_env, ikind, nterms(iatom))
86 22 : CALL get_qs_kind(qs_kind_set(ikind), pao_potentials=pao_potentials)
87 22 : CPASSERT(SIZE(pao_potentials) == 1)
88 54 : maxl = MAX(maxl, pao_potentials(1)%maxl)
89 : END DO
90 10 : CALL init_orbital_pointers(maxl) ! needs to be called before gth_calc_term()
91 :
92 : ! allocate matrix_V_terms
93 10 : CALL dbcsr_get_info(matrix_s(1)%matrix, row_blk_size=blk_sizes_pri)
94 54 : col_blk_size = SUM(nterms)
95 64 : 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 10 : col_blk_size=col_blk_size)
102 10 : CALL dbcsr_reserve_diag_blocks(pao%matrix_V_terms)
103 10 : 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 10 : !$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 10 : IF (pao%precondition) &
126 4 : CALL pao_param_gth_preconditioner(pao, qs_env, nterms)
127 :
128 10 : DEALLOCATE (row_blk_size, col_blk_size, nterms)
129 10 : CALL timestop(handle)
130 10 : END SUBROUTINE pao_param_init_gth
131 :
132 : ! **************************************************************************************************
133 : !> \brief Finalize the GTH potential parametrization
134 : !> \param pao ...
135 : ! **************************************************************************************************
136 10 : SUBROUTINE pao_param_finalize_gth(pao)
137 : TYPE(pao_env_type), POINTER :: pao
138 :
139 10 : CALL dbcsr_release(pao%matrix_V_terms)
140 10 : IF (pao%precondition) THEN
141 4 : CALL dbcsr_release(pao%matrix_precon)
142 4 : CALL dbcsr_release(pao%matrix_precon_inv)
143 : END IF
144 :
145 10 : 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 8 : 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 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_handle)
175 4 : CALL group%set_handle(group_handle)
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)) &
213 3978 : block_overlap = block
214 :
215 120 : 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 4 : threshold=1e-2_dp, converged=arnoldi_converged)
222 6 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| GTH-preconditioner converged, min, max, max/min:", &
223 4 : arnoldi_converged, eval_min, eval_max, eval_max/eval_min
224 :
225 4 : CALL dbcsr_create(pao%matrix_precon, template=matrix_gth_overlap)
226 4 : 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 4 : converged=converged)
234 4 : CALL dbcsr_release(matrix_gth_overlap)
235 :
236 4 : IF (.NOT. converged) &
237 0 : CPABORT("PAO: Sqrt of GTH-preconditioner did not converge.")
238 :
239 4 : CALL timestop(handle)
240 4 : 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 2152 : 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 2152 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
261 : TYPE(dbcsr_type) :: matrix_M, matrix_U
262 :
263 2152 : CALL timeset(routineN, handle)
264 2152 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
265 2152 : CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
266 2152 : CALL dbcsr_reserve_diag_blocks(matrix_U)
267 :
268 : !TODO: move this condition into pao_calc_U, use matrix_N as template
269 2152 : IF (gradient) THEN
270 322 : CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
271 322 : CALL pao_calc_U_gth(pao, matrix_U, matrix_M, pao%matrix_G, penalty)
272 322 : CALL dbcsr_release(matrix_M)
273 : ELSE
274 1830 : CALL pao_calc_U_gth(pao, matrix_U, penalty=penalty)
275 : END IF
276 :
277 2152 : CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
278 2152 : CALL dbcsr_release(matrix_U)
279 2152 : CALL timestop(handle)
280 2152 : 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 2152 : 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 2152 : INTEGER, DIMENSION(:), POINTER :: nterms
301 : LOGICAL :: found
302 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
303 2152 : REAL(dp), DIMENSION(:), POINTER :: world_G, world_X
304 2152 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_M1, block_M2, block_U, &
305 2152 : block_V, block_V_term, block_X, &
306 2152 : vec_V_terms
307 : TYPE(dbcsr_iterator_type) :: iter
308 : TYPE(mp_comm_type) :: group
309 :
310 2152 : CALL timeset(routineN, handle)
311 :
312 2152 : CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nterms, group=group_handle)
313 2152 : CALL group%set_handle(group_handle)
314 2152 : natoms = SIZE(nterms)
315 6456 : ALLOCATE (gaps(natoms))
316 7620 : gaps(:) = HUGE(dp)
317 :
318 : ! allocate arrays for world-view
319 21696 : ALLOCATE (world_X(SUM(nterms)), world_G(SUM(nterms)))
320 202168 : world_X = 0.0_dp; world_G = 0.0_dp
321 :
322 : ! collect world_X from atomic blocks
323 2152 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
324 4886 : DO WHILE (dbcsr_iterator_blocks_left(iter))
325 2734 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
326 2734 : iatom = arow; CPASSERT(arow == acol)
327 4975 : idx = SUM(nterms(1:iatom - 1))
328 104894 : world_X(idx + 1:idx + nterms(iatom)) = block_X(:, 1)
329 : END DO
330 2152 : CALL dbcsr_iterator_stop(iter)
331 202168 : CALL group%sum(world_X) ! sync world view across MPI ranks
332 :
333 : ! loop over atoms
334 2152 : CALL dbcsr_iterator_start(iter, matrix_U)
335 4886 : DO WHILE (dbcsr_iterator_blocks_left(iter))
336 2734 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_U)
337 2734 : iatom = arow; CPASSERT(arow == acol)
338 2734 : n = SIZE(block_U, 1)
339 2734 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=vec_V_terms, found=found)
340 2734 : CPASSERT(ASSOCIATED(vec_V_terms))
341 :
342 : ! calculate potential V of i'th atom
343 10936 : ALLOCATE (block_V(n, n))
344 173370 : block_V = 0.0_dp
345 120226 : DO iterm = 1, SIZE(world_X)
346 117492 : block_V_term(1:n, 1:n) => vec_V_terms(:, iterm) ! map column into matrix
347 12486706 : block_V = block_V + world_X(iterm)*block_V_term
348 : END DO
349 :
350 : ! calculate gradient block of i'th atom
351 2734 : IF (.NOT. PRESENT(matrix_G)) THEN
352 2288 : 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 446 : CPASSERT(PRESENT(matrix_M1))
356 446 : CALL dbcsr_get_block_p(matrix=matrix_M1, row=iatom, col=iatom, block=block_M1, found=found)
357 1338 : ALLOCATE (block_M2(n, n))
358 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
359 446 : M1=block_M1, G=block_M2, gap=gaps(iatom))
360 16910 : DO iterm = 1, SIZE(world_G)
361 16464 : block_V_term(1:n, 1:n) => vec_V_terms(:, iterm) ! map column into matrix
362 1076270 : world_G(iterm) = world_G(iterm) + SUM(block_V_term*block_M2)
363 : END DO
364 892 : DEALLOCATE (block_M2)
365 : END IF
366 10354 : DEALLOCATE (block_V)
367 : END DO
368 2152 : CALL dbcsr_iterator_stop(iter)
369 :
370 : ! distribute world_G across atomic blocks
371 2152 : IF (PRESENT(matrix_G)) THEN
372 25810 : CALL group%sum(world_G) ! sync world view across MPI ranks
373 322 : CALL dbcsr_iterator_start(iter, matrix_G)
374 768 : DO WHILE (dbcsr_iterator_blocks_left(iter))
375 446 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_G)
376 446 : iatom = arow; CPASSERT(arow == acol)
377 855 : idx = SUM(nterms(1:iatom - 1))
378 13512 : block_G(:, 1) = world_G(idx + 1:idx + nterms(iatom))
379 : END DO
380 322 : CALL dbcsr_iterator_stop(iter)
381 : END IF
382 :
383 2152 : DEALLOCATE (world_X, world_G)
384 :
385 : ! sum penalty energies across ranks
386 2152 : IF (PRESENT(penalty)) &
387 2142 : CALL group%sum(penalty)
388 :
389 : ! print homo-lumo gap encountered by fock-layer
390 2152 : CALL group%min(gaps)
391 2152 : IF (pao%iw_gap > 0) THEN
392 2208 : DO iatom = 1, natoms
393 2208 : WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
394 : END DO
395 552 : CALL m_flush(pao%iw_gap)
396 : END IF
397 :
398 : ! one-line summary
399 2152 : IF (pao%iw > 0) THEN
400 8696 : WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", MINVAL(gaps), " for atom:", MINLOC(gaps)
401 : END IF
402 :
403 2152 : DEALLOCATE (gaps)
404 2152 : CALL timestop(handle)
405 :
406 6456 : 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 44 : 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 44 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
421 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
422 :
423 44 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
424 44 : CALL get_qs_kind(qs_kind_set(ikind), pao_potentials=pao_potentials)
425 :
426 44 : IF (SIZE(pao_potentials) /= 1) &
427 0 : CPABORT("GTH parametrization requires exactly one PAO_POTENTIAL section per KIND")
428 :
429 44 : max_projector = pao_potentials(1)%max_projector
430 44 : maxl = pao_potentials(1)%maxl
431 :
432 44 : IF (maxl < 0) &
433 0 : CPABORT("GTH parametrization requires non-negative PAO_POTENTIAL%MAXL")
434 :
435 44 : IF (max_projector < 0) &
436 0 : CPABORT("GTH parametrization requires non-negative PAO_POTENTIAL%MAX_PROJECTOR")
437 :
438 44 : IF (MOD(maxl, 2) /= 0) &
439 0 : CPABORT("GTH parametrization requires even-numbered PAO_POTENTIAL%MAXL")
440 :
441 44 : ncombis = (max_projector + 1)*(max_projector + 2)/2
442 44 : nparams = ncombis*(maxl/2 + 1)
443 :
444 44 : 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 252 : 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 252 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
466 252 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
467 252 : 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 252 : qs_kind_set=qs_kind_set)
473 :
474 : ! get GTH-settings from remote atom
475 252 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
476 252 : CALL get_qs_kind(qs_kind_set(jkind), pao_potentials=pao_potentials)
477 252 : CPASSERT(SIZE(pao_potentials) == 1)
478 252 : pot_max_projector = pao_potentials(1)%max_projector
479 252 : pot_maxl = pao_potentials(1)%maxl
480 252 : pot_beta = pao_potentials(1)%beta
481 :
482 252 : c = 0
483 612 : outer: DO lpot = 0, pot_maxl, 2
484 2252 : DO max_l = 0, pot_max_projector
485 4718 : DO min_l = 0, max_l
486 2970 : c = c + 1
487 4358 : IF (c == kterm) EXIT outer
488 : END DO
489 : END DO
490 : END DO outer
491 :
492 : ! get basis-set of central atom
493 252 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
494 252 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set)
495 :
496 1008 : Ra = particle_set(iatom)%r
497 1008 : Rb = particle_set(jatom)%r
498 252 : Rab = pbc(ra, rb, cell)
499 :
500 15108 : block_V = 0.0_dp
501 : CALL pao_calc_gaussian(basis_set, block_V, Rab=Rab, lpot=lpot, &
502 252 : min_l=min_l, max_l=max_l, beta=pot_beta, weight=1.0_dp)
503 :
504 252 : END SUBROUTINE gth_calc_term
505 :
506 : ! **************************************************************************************************
507 : !> \brief Calculate initial guess for matrix_X
508 : !> \param pao ...
509 : ! **************************************************************************************************
510 10 : SUBROUTINE pao_param_initguess_gth(pao)
511 : TYPE(pao_env_type), POINTER :: pao
512 :
513 : INTEGER :: acol, arow
514 10 : REAL(dp), DIMENSION(:, :), POINTER :: block_X
515 : TYPE(dbcsr_iterator_type) :: iter
516 :
517 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
518 10 : !$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 10 : END SUBROUTINE pao_param_initguess_gth
533 :
534 : END MODULE pao_param_gth
|
