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