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 Common framework for a linear parametrization of the potential.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_param_linpot
13 : USE atomic_kind_types, ONLY: get_atomic_kind
14 : USE basis_set_types, ONLY: gto_basis_set_type
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: &
17 : dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
18 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
19 : dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type
20 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
21 : USE kinds, ONLY: dp
22 : USE machine, ONLY: m_flush
23 : USE mathlib, ONLY: diamat_all
24 : USE message_passing, ONLY: mp_comm_type,&
25 : mp_para_env_type
26 : USE pao_input, ONLY: pao_fock_param,&
27 : pao_rotinv_param
28 : USE pao_linpot_full, ONLY: linpot_full_calc_terms,&
29 : linpot_full_count_terms
30 : USE pao_linpot_rotinv, ONLY: linpot_rotinv_calc_forces,&
31 : linpot_rotinv_calc_terms,&
32 : linpot_rotinv_count_terms
33 : USE pao_param_fock, ONLY: pao_calc_U_block_fock
34 : USE pao_param_methods, ONLY: pao_calc_AB_from_U,&
35 : pao_calc_grad_lnv_wrt_U
36 : USE pao_potentials, ONLY: pao_guess_initial_potential
37 : USE pao_types, ONLY: pao_env_type
38 : USE particle_types, ONLY: particle_type
39 : USE qs_environment_types, ONLY: get_qs_env,&
40 : qs_environment_type
41 : USE qs_kind_types, ONLY: get_qs_kind,&
42 : qs_kind_type
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : PRIVATE
48 :
49 : PUBLIC :: pao_param_init_linpot, pao_param_finalize_linpot, pao_calc_AB_linpot
50 : PUBLIC :: pao_param_count_linpot, pao_param_initguess_linpot
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief Initialize the linear potential parametrization
56 : !> \param pao ...
57 : !> \param qs_env ...
58 : ! **************************************************************************************************
59 234 : SUBROUTINE pao_param_init_linpot(pao, qs_env)
60 : TYPE(pao_env_type), POINTER :: pao
61 : TYPE(qs_environment_type), POINTER :: qs_env
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_init_linpot'
64 :
65 : INTEGER :: acol, arow, handle, iatom, ikind, N, &
66 : natoms, nterms
67 234 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pri, col_blk_size, row_blk_size
68 234 : REAL(dp), DIMENSION(:, :), POINTER :: block_V_terms
69 234 : REAL(dp), DIMENSION(:, :, :), POINTER :: V_blocks
70 : TYPE(dbcsr_iterator_type) :: iter
71 : TYPE(dft_control_type), POINTER :: dft_control
72 : TYPE(mp_para_env_type), POINTER :: para_env
73 234 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
74 :
75 234 : CALL timeset(routineN, handle)
76 :
77 : CALL get_qs_env(qs_env, &
78 : para_env=para_env, &
79 : dft_control=dft_control, &
80 : particle_set=particle_set, &
81 234 : natom=natoms)
82 :
83 234 : IF (dft_control%nspins /= 1) CPABORT("open shell not yet implemented")
84 :
85 : ! figure out number of potential terms
86 936 : ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
87 714 : DO iatom = 1, natoms
88 480 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
89 480 : CALL pao_param_count_linpot(pao, qs_env, ikind, nterms)
90 714 : col_blk_size(iatom) = nterms
91 : END DO
92 :
93 : ! allocate matrix_V_terms
94 234 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri)
95 1428 : 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 234 : col_blk_size=col_blk_size)
102 234 : CALL dbcsr_reserve_diag_blocks(pao%matrix_V_terms)
103 234 : DEALLOCATE (row_blk_size, col_blk_size)
104 :
105 : ! calculate, normalize, and store potential terms as rows of block_V_terms
106 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,blk_sizes_pri) &
107 234 : !$OMP PRIVATE(iter,arow,acol,iatom,N,nterms,block_V_terms,V_blocks)
108 : CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
109 : DO WHILE (dbcsr_iterator_blocks_left(iter))
110 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_V_terms)
111 : iatom = arow; CPASSERT(arow == acol)
112 : nterms = SIZE(block_V_terms, 2)
113 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
114 : N = blk_sizes_pri(iatom)
115 : CPASSERT(N*N == SIZE(block_V_terms, 1))
116 : ALLOCATE (V_blocks(N, N, nterms))
117 : CALL linpot_calc_terms(pao, qs_env, iatom, V_blocks)
118 : block_V_terms = RESHAPE(V_blocks, (/N*N, nterms/)) ! convert matrices into vectors
119 : DEALLOCATE (V_blocks)
120 : END DO
121 : CALL dbcsr_iterator_stop(iter)
122 : !$OMP END PARALLEL
123 :
124 234 : CALL pao_param_linpot_regularizer(pao)
125 :
126 234 : IF (pao%precondition) &
127 12 : CALL pao_param_linpot_preconditioner(pao)
128 :
129 234 : CALL para_env%sync() ! ensure that timestop is not called too early
130 :
131 234 : CALL timestop(handle)
132 234 : END SUBROUTINE pao_param_init_linpot
133 :
134 : ! **************************************************************************************************
135 : !> \brief Builds the regularization metric matrix_R
136 : !> \param pao ...
137 : ! **************************************************************************************************
138 234 : SUBROUTINE pao_param_linpot_regularizer(pao)
139 : TYPE(pao_env_type), POINTER :: pao
140 :
141 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_linpot_regularizer'
142 :
143 : INTEGER :: acol, arow, handle, i, iatom, j, k, &
144 : nterms
145 234 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
146 : LOGICAL :: found
147 : REAL(dp) :: v, w
148 234 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
149 234 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs
150 234 : REAL(dp), DIMENSION(:, :), POINTER :: block_R, V_terms
151 : TYPE(dbcsr_iterator_type) :: iter
152 :
153 234 : CALL timeset(routineN, handle)
154 :
155 234 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot regularizer"
156 :
157 234 : CALL dbcsr_get_info(pao%matrix_V_terms, col_blk_size=blk_sizes_nterms)
158 :
159 : ! build regularization metric
160 : CALL dbcsr_create(pao%matrix_R, &
161 : template=pao%matrix_V_terms, &
162 : matrix_type="N", &
163 : row_blk_size=blk_sizes_nterms, &
164 : col_blk_size=blk_sizes_nterms, &
165 234 : name="PAO matrix_R")
166 234 : CALL dbcsr_reserve_diag_blocks(pao%matrix_R)
167 :
168 : ! fill matrix_R
169 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
170 234 : !$OMP PRIVATE(iter,arow,acol,iatom,block_R,V_terms,found,nterms,S,S_evecs,S_evals,k,i,j,v,w)
171 : CALL dbcsr_iterator_start(iter, pao%matrix_R)
172 : DO WHILE (dbcsr_iterator_blocks_left(iter))
173 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_R)
174 : iatom = arow; CPASSERT(arow == acol)
175 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=V_terms, found=found)
176 : CPASSERT(ASSOCIATED(V_terms))
177 : nterms = SIZE(V_terms, 2)
178 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
179 :
180 : ! build overlap matrix
181 : ALLOCATE (S(nterms, nterms))
182 : S(:, :) = MATMUL(TRANSPOSE(V_terms), V_terms)
183 :
184 : ! diagonalize S
185 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
186 : S_evecs(:, :) = S
187 : CALL diamat_all(S_evecs, S_evals)
188 :
189 : block_R = 0.0_dp
190 : DO k = 1, nterms
191 : v = pao%linpot_regu_delta/S_evals(k)
192 : w = pao%linpot_regu_strength*MIN(1.0_dp, ABS(v))
193 : DO i = 1, nterms
194 : DO j = 1, nterms
195 : block_R(i, j) = block_R(i, j) + w*S_evecs(i, k)*S_evecs(j, k)
196 : END DO
197 : END DO
198 : END DO
199 :
200 : ! clean up
201 : DEALLOCATE (S, S_evals, S_evecs)
202 : END DO
203 : CALL dbcsr_iterator_stop(iter)
204 : !$OMP END PARALLEL
205 :
206 234 : CALL timestop(handle)
207 468 : END SUBROUTINE pao_param_linpot_regularizer
208 :
209 : ! **************************************************************************************************
210 : !> \brief Builds the preconditioner matrix_precon and matrix_precon_inv
211 : !> \param pao ...
212 : ! **************************************************************************************************
213 12 : SUBROUTINE pao_param_linpot_preconditioner(pao)
214 : TYPE(pao_env_type), POINTER :: pao
215 :
216 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_linpot_preconditioner'
217 :
218 : INTEGER :: acol, arow, handle, i, iatom, j, k, &
219 : nterms
220 12 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_nterms
221 : LOGICAL :: found
222 : REAL(dp) :: eval_capped
223 12 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
224 12 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs
225 12 : REAL(dp), DIMENSION(:, :), POINTER :: block_precon, block_precon_inv, &
226 12 : block_V_terms
227 : TYPE(dbcsr_iterator_type) :: iter
228 :
229 12 : CALL timeset(routineN, handle)
230 :
231 12 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Building linpot preconditioner"
232 :
233 12 : CALL dbcsr_get_info(pao%matrix_V_terms, col_blk_size=blk_sizes_nterms)
234 :
235 : CALL dbcsr_create(pao%matrix_precon, &
236 : template=pao%matrix_V_terms, &
237 : matrix_type="N", &
238 : row_blk_size=blk_sizes_nterms, &
239 : col_blk_size=blk_sizes_nterms, &
240 12 : name="PAO matrix_precon")
241 12 : CALL dbcsr_reserve_diag_blocks(pao%matrix_precon)
242 :
243 12 : CALL dbcsr_create(pao%matrix_precon_inv, template=pao%matrix_precon, name="PAO matrix_precon_inv")
244 12 : CALL dbcsr_reserve_diag_blocks(pao%matrix_precon_inv)
245 :
246 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
247 12 : !$OMP PRIVATE(iter,arow,acol,iatom,block_V_terms,block_precon,block_precon_inv,found,nterms,S,S_evals,S_evecs,i,j,k,eval_capped)
248 : CALL dbcsr_iterator_start(iter, pao%matrix_V_terms)
249 : DO WHILE (dbcsr_iterator_blocks_left(iter))
250 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_V_terms)
251 : iatom = arow; CPASSERT(arow == acol)
252 : nterms = SIZE(block_V_terms, 2)
253 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
254 :
255 : CALL dbcsr_get_block_p(matrix=pao%matrix_precon, row=iatom, col=iatom, block=block_precon, found=found)
256 : CALL dbcsr_get_block_p(matrix=pao%matrix_precon_inv, row=iatom, col=iatom, block=block_precon_inv, found=found)
257 : CPASSERT(ASSOCIATED(block_precon))
258 : CPASSERT(ASSOCIATED(block_precon_inv))
259 :
260 : ALLOCATE (S(nterms, nterms))
261 : S(:, :) = MATMUL(TRANSPOSE(block_V_terms), block_V_terms)
262 :
263 : ! diagonalize S
264 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
265 : S_evecs(:, :) = S
266 : CALL diamat_all(S_evecs, S_evals)
267 :
268 : ! construct 1/Sqrt(S) and Sqrt(S)
269 : block_precon = 0.0_dp
270 : block_precon_inv = 0.0_dp
271 : DO k = 1, nterms
272 : eval_capped = MAX(pao%linpot_precon_delta, S_evals(k)) ! too small eigenvalues are hurtful
273 : DO i = 1, nterms
274 : DO j = 1, nterms
275 : block_precon(i, j) = block_precon(i, j) + S_evecs(i, k)*S_evecs(j, k)/SQRT(eval_capped)
276 : block_precon_inv(i, j) = block_precon_inv(i, j) + S_evecs(i, k)*S_evecs(j, k)*SQRT(eval_capped)
277 : END DO
278 : END DO
279 : END DO
280 :
281 : DEALLOCATE (S, S_evecs, S_evals)
282 : END DO
283 : CALL dbcsr_iterator_stop(iter)
284 : !$OMP END PARALLEL
285 :
286 12 : CALL timestop(handle)
287 24 : END SUBROUTINE pao_param_linpot_preconditioner
288 :
289 : ! **************************************************************************************************
290 : !> \brief Finalize the linear potential parametrization
291 : !> \param pao ...
292 : ! **************************************************************************************************
293 234 : SUBROUTINE pao_param_finalize_linpot(pao)
294 : TYPE(pao_env_type), POINTER :: pao
295 :
296 234 : CALL dbcsr_release(pao%matrix_V_terms)
297 234 : CALL dbcsr_release(pao%matrix_R)
298 :
299 234 : IF (pao%precondition) THEN
300 12 : CALL dbcsr_release(pao%matrix_precon)
301 12 : CALL dbcsr_release(pao%matrix_precon_inv)
302 : END IF
303 :
304 234 : END SUBROUTINE pao_param_finalize_linpot
305 :
306 : ! **************************************************************************************************
307 : !> \brief Returns the number of potential terms for given atomic kind
308 : !> \param pao ...
309 : !> \param qs_env ...
310 : !> \param ikind ...
311 : !> \param nparams ...
312 : ! **************************************************************************************************
313 1344 : SUBROUTINE pao_param_count_linpot(pao, qs_env, ikind, nparams)
314 : TYPE(pao_env_type), POINTER :: pao
315 : TYPE(qs_environment_type), POINTER :: qs_env
316 : INTEGER, INTENT(IN) :: ikind
317 : INTEGER, INTENT(OUT) :: nparams
318 :
319 : INTEGER :: pao_basis_size
320 : TYPE(gto_basis_set_type), POINTER :: basis_set
321 672 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
322 :
323 672 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
324 :
325 : CALL get_qs_kind(qs_kind_set(ikind), &
326 : basis_set=basis_set, &
327 672 : pao_basis_size=pao_basis_size)
328 :
329 672 : IF (pao_basis_size == basis_set%nsgf) THEN
330 26 : nparams = 0 ! pao disabled for iatom
331 :
332 : ELSE
333 754 : SELECT CASE (pao%parameterization)
334 : CASE (pao_fock_param)
335 646 : CALL linpot_full_count_terms(qs_env, ikind, nterms=nparams)
336 : CASE (pao_rotinv_param)
337 538 : CALL linpot_rotinv_count_terms(qs_env, ikind, nterms=nparams)
338 : CASE DEFAULT
339 646 : CPABORT("unknown parameterization")
340 : END SELECT
341 : END IF
342 :
343 672 : END SUBROUTINE pao_param_count_linpot
344 :
345 : ! **************************************************************************************************
346 : !> \brief Takes current matrix_X and calculates the matrices A and B.
347 : !> \param pao ...
348 : !> \param qs_env ...
349 : !> \param ls_scf_env ...
350 : !> \param gradient ...
351 : !> \param penalty ...
352 : !> \param forces ...
353 : ! **************************************************************************************************
354 8192 : SUBROUTINE pao_calc_AB_linpot(pao, qs_env, ls_scf_env, gradient, penalty, forces)
355 : TYPE(pao_env_type), POINTER :: pao
356 : TYPE(qs_environment_type), POINTER :: qs_env
357 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
358 : LOGICAL, INTENT(IN) :: gradient
359 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
360 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
361 :
362 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_AB_linpot'
363 :
364 : INTEGER :: handle
365 8192 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
366 : TYPE(dbcsr_type) :: matrix_M, matrix_U
367 :
368 8192 : CALL timeset(routineN, handle)
369 8192 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
370 8192 : CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
371 8192 : CALL dbcsr_reserve_diag_blocks(matrix_U)
372 :
373 : !TODO: move this condition into pao_calc_U, use matrix_N as template
374 8192 : IF (gradient) THEN
375 1616 : CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
376 3198 : CALL pao_calc_U_linpot(pao, qs_env, matrix_U, matrix_M, pao%matrix_G, penalty, forces)
377 1616 : CALL dbcsr_release(matrix_M)
378 : ELSE
379 6576 : CALL pao_calc_U_linpot(pao, qs_env, matrix_U, penalty=penalty)
380 : END IF
381 :
382 8192 : CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
383 8192 : CALL dbcsr_release(matrix_U)
384 8192 : CALL timestop(handle)
385 8192 : END SUBROUTINE pao_calc_AB_linpot
386 :
387 : ! **************************************************************************************************
388 : !> \brief Calculate new matrix U and optinally its gradient G
389 : !> \param pao ...
390 : !> \param qs_env ...
391 : !> \param matrix_U ...
392 : !> \param matrix_M ...
393 : !> \param matrix_G ...
394 : !> \param penalty ...
395 : !> \param forces ...
396 : ! **************************************************************************************************
397 8192 : SUBROUTINE pao_calc_U_linpot(pao, qs_env, matrix_U, matrix_M, matrix_G, penalty, forces)
398 : TYPE(pao_env_type), POINTER :: pao
399 : TYPE(qs_environment_type), POINTER :: qs_env
400 : TYPE(dbcsr_type) :: matrix_U
401 : TYPE(dbcsr_type), OPTIONAL :: matrix_M, matrix_G
402 : REAL(dp), INTENT(INOUT), OPTIONAL :: penalty
403 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
404 :
405 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_U_linpot'
406 :
407 : INTEGER :: acol, arow, handle, iatom, kterm, n, &
408 : natoms, nterms
409 : LOGICAL :: found
410 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
411 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: evals
412 8192 : REAL(dp), DIMENSION(:), POINTER :: vec_M2, vec_V
413 8192 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_M1, block_M2, block_R, &
414 8192 : block_U, block_V, block_V_terms, &
415 8192 : block_X
416 8192 : REAL(dp), DIMENSION(:, :, :), POINTER :: M_blocks
417 : REAL(KIND=dp) :: regu_energy
418 : TYPE(dbcsr_iterator_type) :: iter
419 : TYPE(mp_comm_type) :: group
420 :
421 8192 : CALL timeset(routineN, handle)
422 :
423 8192 : CPASSERT(PRESENT(matrix_G) .EQV. PRESENT(matrix_M))
424 :
425 8192 : CALL get_qs_env(qs_env, natom=natoms)
426 40960 : ALLOCATE (gaps(natoms), evals(10, natoms)) ! printing 10 eigenvalues seems reasonable
427 244604 : evals(:, :) = 0.0_dp
428 29684 : gaps(:) = HUGE(1.0_dp)
429 8192 : regu_energy = 0.0_dp
430 8192 : CALL dbcsr_get_info(matrix_U, group=group)
431 :
432 8192 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
433 18938 : DO WHILE (dbcsr_iterator_blocks_left(iter))
434 10746 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
435 10746 : iatom = arow; CPASSERT(arow == acol)
436 10746 : CALL dbcsr_get_block_p(matrix=pao%matrix_R, row=iatom, col=iatom, block=block_R, found=found)
437 10746 : CALL dbcsr_get_block_p(matrix=matrix_U, row=iatom, col=iatom, block=block_U, found=found)
438 10746 : CPASSERT(ASSOCIATED(block_R) .AND. ASSOCIATED(block_U))
439 10746 : n = SIZE(block_U, 1)
440 :
441 : ! calculate potential V
442 32238 : ALLOCATE (vec_V(n*n))
443 647165 : vec_V(:) = 0.0_dp
444 10746 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=block_V_terms, found=found)
445 10746 : CPASSERT(ASSOCIATED(block_V_terms))
446 10746 : nterms = SIZE(block_V_terms, 2)
447 10746 : IF (nterms > 0) & ! protect against corner-case of zero pao parameters
448 71735320 : vec_V = MATMUL(block_V_terms, block_X(:, 1))
449 10746 : block_V(1:n, 1:n) => vec_V(:) ! map vector into matrix
450 :
451 : ! symmetrize
452 1431894 : IF (MAXVAL(ABS(block_V - TRANSPOSE(block_V))/MAX(1.0_dp, MAXVAL(ABS(block_V)))) > 1e-12) &
453 0 : CPABORT("block_V not symmetric")
454 1442640 : block_V = 0.5_dp*(block_V + TRANSPOSE(block_V)) ! symmetrize exactly
455 :
456 : ! regularization energy
457 : ! protect against corner-case of zero pao parameters
458 10746 : IF (PRESENT(penalty) .AND. nterms > 0) &
459 30347076 : regu_energy = regu_energy + DOT_PRODUCT(block_X(:, 1), MATMUL(block_R, block_X(:, 1)))
460 :
461 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
462 10746 : gap=gaps(iatom), evals=evals(:, iatom))
463 :
464 10746 : IF (PRESENT(matrix_G)) THEN ! TURNING POINT (if calc grad) --------------------------------
465 2132 : CPASSERT(PRESENT(matrix_M))
466 2132 : CALL dbcsr_get_block_p(matrix=matrix_M, row=iatom, col=iatom, block=block_M1, found=found)
467 :
468 : ! corner-cases: block_M1 might have been filtered out or there might be zero pao parameters
469 6396 : IF (ASSOCIATED(block_M1) .AND. SIZE(block_V_terms) > 0) THEN
470 4038 : ALLOCATE (vec_M2(n*n))
471 2019 : block_M2(1:n, 1:n) => vec_M2(:) ! map vector into matrix
472 : !TODO: this 2nd call does double work. However, *sometimes* this branch is not taken.
473 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
474 2019 : M1=block_M1, G=block_M2, gap=gaps(iatom), evals=evals(:, iatom))
475 124173 : IF (MAXVAL(ABS(block_M2 - TRANSPOSE(block_M2))) > 1e-14_dp) &
476 0 : CPABORT("matrix not symmetric")
477 :
478 : ! gradient dE/dX
479 2019 : IF (PRESENT(matrix_G)) THEN
480 2019 : CALL dbcsr_get_block_p(matrix=matrix_G, row=iatom, col=iatom, block=block_G, found=found)
481 2019 : CPASSERT(ASSOCIATED(block_G))
482 6582106 : block_G(:, 1) = MATMUL(vec_M2, block_V_terms)
483 2019 : IF (PRESENT(penalty)) &
484 10731904 : block_G = block_G + 2.0_dp*MATMUL(block_R, block_X) ! regularization gradient
485 : END IF
486 :
487 : ! forced dE/dR
488 2019 : IF (PRESENT(forces)) THEN
489 170 : ALLOCATE (M_blocks(n, n, nterms))
490 296 : DO kterm = 1, nterms
491 16806 : M_blocks(:, :, kterm) = block_M2*block_X(kterm, 1)
492 : END DO
493 34 : CALL linpot_calc_forces(pao, qs_env, iatom=iatom, M_blocks=M_blocks, forces=forces)
494 34 : DEALLOCATE (M_blocks)
495 : END IF
496 :
497 2019 : DEALLOCATE (vec_M2)
498 : END IF
499 : END IF
500 51176 : DEALLOCATE (vec_V)
501 : END DO
502 8192 : CALL dbcsr_iterator_stop(iter)
503 :
504 8192 : IF (PRESENT(penalty)) THEN
505 : ! sum penalty energies across ranks
506 7924 : CALL group%sum(penalty)
507 7924 : CALL group%sum(regu_energy)
508 7924 : penalty = penalty + regu_energy
509 : END IF
510 :
511 : ! print stuff, but not during second invocation for forces
512 8192 : IF (.NOT. PRESENT(forces)) THEN
513 : ! print eigenvalues from fock-layer
514 8158 : CALL group%sum(evals)
515 8158 : IF (pao%iw_fockev > 0) THEN
516 2000 : DO iatom = 1, natoms
517 2000 : WRITE (pao%iw_fockev, *) "PAO| atom:", iatom, " fock evals around gap:", evals(:, iatom)
518 : END DO
519 500 : CALL m_flush(pao%iw_fockev)
520 : END IF
521 : ! print homo-lumo gap encountered by fock-layer
522 8158 : CALL group%min(gaps)
523 8158 : IF (pao%iw_gap > 0) THEN
524 2000 : DO iatom = 1, natoms
525 2000 : WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
526 : END DO
527 : END IF
528 : ! one-line summaries
529 8158 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| linpot regularization energy:", regu_energy
530 37738 : IF (pao%iw > 0) WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", MINVAL(gaps), " for atom:", MINLOC(gaps)
531 : END IF
532 :
533 8192 : DEALLOCATE (gaps, evals)
534 8192 : CALL timestop(handle)
535 :
536 16384 : END SUBROUTINE pao_calc_U_linpot
537 :
538 : ! **************************************************************************************************
539 : !> \brief Internal routine, calculates terms in potential parametrization
540 : !> \param pao ...
541 : !> \param qs_env ...
542 : !> \param iatom ...
543 : !> \param V_blocks ...
544 : ! **************************************************************************************************
545 234 : SUBROUTINE linpot_calc_terms(pao, qs_env, iatom, V_blocks)
546 : TYPE(pao_env_type), POINTER :: pao
547 : TYPE(qs_environment_type), POINTER :: qs_env
548 : INTEGER, INTENT(IN) :: iatom
549 : REAL(dp), DIMENSION(:, :, :), INTENT(OUT) :: V_blocks
550 :
551 273 : SELECT CASE (pao%parameterization)
552 : CASE (pao_fock_param)
553 39 : CALL linpot_full_calc_terms(V_blocks)
554 : CASE (pao_rotinv_param)
555 195 : CALL linpot_rotinv_calc_terms(qs_env, iatom, V_blocks)
556 : CASE DEFAULT
557 234 : CPABORT("unknown parameterization")
558 : END SELECT
559 :
560 234 : END SUBROUTINE linpot_calc_terms
561 :
562 : ! **************************************************************************************************
563 : !> \brief Internal routine, calculates force contributions from potential parametrization
564 : !> \param pao ...
565 : !> \param qs_env ...
566 : !> \param iatom ...
567 : !> \param M_blocks ...
568 : !> \param forces ...
569 : ! **************************************************************************************************
570 34 : SUBROUTINE linpot_calc_forces(pao, qs_env, iatom, M_blocks, forces)
571 : TYPE(pao_env_type), POINTER :: pao
572 : TYPE(qs_environment_type), POINTER :: qs_env
573 : INTEGER, INTENT(IN) :: iatom
574 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: M_blocks
575 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: forces
576 :
577 66 : SELECT CASE (pao%parameterization)
578 : CASE (pao_fock_param)
579 : ! no force contributions
580 : CASE (pao_rotinv_param)
581 32 : CALL linpot_rotinv_calc_forces(qs_env, iatom, M_blocks, forces)
582 : CASE DEFAULT
583 34 : CPABORT("unknown parameterization")
584 : END SELECT
585 :
586 34 : END SUBROUTINE linpot_calc_forces
587 :
588 : ! **************************************************************************************************
589 : !> \brief Calculate initial guess for matrix_X
590 : !> \param pao ...
591 : !> \param qs_env ...
592 : ! **************************************************************************************************
593 34 : SUBROUTINE pao_param_initguess_linpot(pao, qs_env)
594 : TYPE(pao_env_type), POINTER :: pao
595 : TYPE(qs_environment_type), POINTER :: qs_env
596 :
597 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_initguess_linpot'
598 :
599 : INTEGER :: acol, arow, handle, i, iatom, j, k, n, &
600 : nterms
601 34 : INTEGER, DIMENSION(:), POINTER :: pri_basis_size
602 : LOGICAL :: found
603 : REAL(dp) :: w
604 34 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
605 34 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs, S_inv
606 34 : REAL(dp), DIMENSION(:), POINTER :: V_guess_vec
607 34 : REAL(dp), DIMENSION(:, :), POINTER :: block_X, V_guess, V_terms
608 : TYPE(dbcsr_iterator_type) :: iter
609 :
610 34 : CALL timeset(routineN, handle)
611 :
612 34 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis_size)
613 :
614 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,pri_basis_size) &
615 34 : !$OMP PRIVATE(iter,arow,acol,iatom,block_X,N,nterms,V_terms,found,V_guess,V_guess_vec,S,S_evecs,S_evals,S_inv,k,i,j,w)
616 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
617 : DO WHILE (dbcsr_iterator_blocks_left(iter))
618 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
619 : iatom = arow; CPASSERT(arow == acol)
620 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=V_terms, found=found)
621 : CPASSERT(ASSOCIATED(V_terms))
622 : nterms = SIZE(V_terms, 2)
623 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
624 :
625 : ! guess initial potential
626 : N = pri_basis_size(iatom)
627 : ALLOCATE (V_guess_vec(n*n))
628 : V_guess(1:n, 1:n) => V_guess_vec
629 : CALL pao_guess_initial_potential(qs_env, iatom, V_guess)
630 :
631 : ! build overlap matrix
632 : ALLOCATE (S(nterms, nterms))
633 : S(:, :) = MATMUL(TRANSPOSE(V_terms), V_terms)
634 :
635 : ! diagonalize S
636 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
637 : S_evecs(:, :) = S
638 : CALL diamat_all(S_evecs, S_evals)
639 :
640 : ! calculate Tikhonov regularized inverse
641 : ALLOCATE (S_inv(nterms, nterms))
642 : S_inv(:, :) = 0.0_dp
643 : DO k = 1, nterms
644 : w = S_evals(k)/(S_evals(k)**2 + pao%linpot_init_delta)
645 : DO i = 1, nterms
646 : DO j = 1, nterms
647 : S_inv(i, j) = S_inv(i, j) + w*S_evecs(i, k)*S_evecs(j, k)
648 : END DO
649 : END DO
650 : END DO
651 :
652 : ! perform fit
653 : block_X(:, 1) = MATMUL(MATMUL(S_inv, TRANSPOSE(V_terms)), V_guess_vec)
654 :
655 : ! clean up
656 : DEALLOCATE (V_guess_vec, S, S_evecs, S_evals, S_inv)
657 : END DO
658 : CALL dbcsr_iterator_stop(iter)
659 : !$OMP END PARALLEL
660 :
661 34 : CALL timestop(handle)
662 68 : END SUBROUTINE pao_param_initguess_linpot
663 :
664 24204 : END MODULE pao_param_linpot
|
