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 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 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 1170 : 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 108 : 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 8210 : 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 8210 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
366 : TYPE(dbcsr_type) :: matrix_M, matrix_U
367 :
368 8210 : CALL timeset(routineN, handle)
369 8210 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
370 8210 : CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
371 8210 : CALL dbcsr_reserve_diag_blocks(matrix_U)
372 :
373 : !TODO: move this condition into pao_calc_U, use matrix_N as template
374 8210 : IF (gradient) THEN
375 1622 : CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
376 3210 : CALL pao_calc_U_linpot(pao, qs_env, matrix_U, matrix_M, pao%matrix_G, penalty, forces)
377 1622 : CALL dbcsr_release(matrix_M)
378 : ELSE
379 6588 : CALL pao_calc_U_linpot(pao, qs_env, matrix_U, penalty=penalty)
380 : END IF
381 :
382 8210 : CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
383 8210 : CALL dbcsr_release(matrix_U)
384 8210 : CALL timestop(handle)
385 8210 : 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 8210 : 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, group_handle, handle, iatom, &
408 : kterm, n, natoms, nterms
409 : LOGICAL :: found
410 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: gaps
411 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: evals
412 8210 : REAL(dp), DIMENSION(:), POINTER :: vec_M2, vec_V
413 8210 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_M1, block_M2, block_R, &
414 8210 : block_U, block_V, block_V_terms, &
415 8210 : block_X
416 8210 : 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 8210 : CALL timeset(routineN, handle)
422 :
423 8210 : CPASSERT(PRESENT(matrix_G) .EQV. PRESENT(matrix_M))
424 :
425 8210 : CALL get_qs_env(qs_env, natom=natoms)
426 41050 : ALLOCATE (gaps(natoms), evals(10, natoms)) ! printing 10 eigenvalues seems reasonable
427 244952 : evals(:, :) = 0.0_dp
428 29732 : gaps(:) = HUGE(1.0_dp)
429 8210 : regu_energy = 0.0_dp
430 8210 : CALL dbcsr_get_info(matrix_U, group=group_handle)
431 8210 : CALL group%set_handle(group_handle)
432 :
433 8210 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
434 18971 : DO WHILE (dbcsr_iterator_blocks_left(iter))
435 10761 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
436 10761 : iatom = arow; CPASSERT(arow == acol)
437 10761 : CALL dbcsr_get_block_p(matrix=pao%matrix_R, row=iatom, col=iatom, block=block_R, found=found)
438 10761 : CALL dbcsr_get_block_p(matrix=matrix_U, row=iatom, col=iatom, block=block_U, found=found)
439 10761 : CPASSERT(ASSOCIATED(block_R) .AND. ASSOCIATED(block_U))
440 10761 : n = SIZE(block_U, 1)
441 :
442 : ! calculate potential V
443 32283 : ALLOCATE (vec_V(n*n))
444 647123 : vec_V(:) = 0.0_dp
445 10761 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=block_V_terms, found=found)
446 10761 : CPASSERT(ASSOCIATED(block_V_terms))
447 10761 : nterms = SIZE(block_V_terms, 2)
448 10761 : IF (nterms > 0) & ! protect against corner-case of zero pao parameters
449 71703382 : vec_V = MATMUL(block_V_terms, block_X(:, 1))
450 10761 : block_V(1:n, 1:n) => vec_V(:) ! map vector into matrix
451 :
452 : ! symmetrize
453 1431897 : IF (MAXVAL(ABS(block_V - TRANSPOSE(block_V))/MAX(1.0_dp, MAXVAL(ABS(block_V)))) > 1e-12) &
454 0 : CPABORT("block_V not symmetric")
455 1431897 : block_V = 0.5_dp*(block_V + TRANSPOSE(block_V)) ! symmetrize exactly
456 :
457 : ! regularization energy
458 : ! protect against corner-case of zero pao parameters
459 10761 : IF (PRESENT(penalty) .AND. nterms > 0) &
460 30339810 : regu_energy = regu_energy + DOT_PRODUCT(block_X(:, 1), MATMUL(block_R, block_X(:, 1)))
461 :
462 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
463 10761 : gap=gaps(iatom), evals=evals(:, iatom))
464 :
465 10761 : IF (PRESENT(matrix_G)) THEN ! TURNING POINT (if calc grad) --------------------------------
466 2138 : CPASSERT(PRESENT(matrix_M))
467 2138 : CALL dbcsr_get_block_p(matrix=matrix_M, row=iatom, col=iatom, block=block_M1, found=found)
468 :
469 : ! corner-cases: block_M1 might have been filtered out or there might be zero pao parameters
470 2138 : IF (ASSOCIATED(block_M1) .AND. SIZE(block_V_terms) > 0) THEN
471 6075 : ALLOCATE (vec_M2(n*n))
472 2025 : block_M2(1:n, 1:n) => vec_M2(:) ! map vector into matrix
473 : !TODO: this 2nd call does double work. However, *sometimes* this branch is not taken.
474 : CALL pao_calc_U_block_fock(pao, iatom=iatom, penalty=penalty, V=block_V, U=block_U, &
475 2025 : M1=block_M1, G=block_M2, gap=gaps(iatom), evals=evals(:, iatom))
476 124359 : IF (MAXVAL(ABS(block_M2 - TRANSPOSE(block_M2))) > 1e-14_dp) &
477 0 : CPABORT("matrix not symmetric")
478 :
479 : ! gradient dE/dX
480 2025 : IF (PRESENT(matrix_G)) THEN
481 2025 : CALL dbcsr_get_block_p(matrix=matrix_G, row=iatom, col=iatom, block=block_G, found=found)
482 2025 : CPASSERT(ASSOCIATED(block_G))
483 6592174 : block_G(:, 1) = MATMUL(vec_M2, block_V_terms)
484 2025 : IF (PRESENT(penalty)) &
485 10738337 : block_G = block_G + 2.0_dp*MATMUL(block_R, block_X) ! regularization gradient
486 : END IF
487 :
488 : ! forced dE/dR
489 2025 : IF (PRESENT(forces)) THEN
490 170 : ALLOCATE (M_blocks(n, n, nterms))
491 296 : DO kterm = 1, nterms
492 16544 : M_blocks(:, :, kterm) = block_M2*block_X(kterm, 1)
493 : END DO
494 34 : CALL linpot_calc_forces(pao, qs_env, iatom=iatom, M_blocks=M_blocks, forces=forces)
495 34 : DEALLOCATE (M_blocks)
496 : END IF
497 :
498 2025 : DEALLOCATE (vec_M2)
499 : END IF
500 : END IF
501 51254 : DEALLOCATE (vec_V)
502 : END DO
503 8210 : CALL dbcsr_iterator_stop(iter)
504 :
505 8210 : IF (PRESENT(penalty)) THEN
506 : ! sum penalty energies across ranks
507 7942 : CALL group%sum(penalty)
508 7942 : CALL group%sum(regu_energy)
509 7942 : penalty = penalty + regu_energy
510 : END IF
511 :
512 : ! print stuff, but not during second invocation for forces
513 8210 : IF (.NOT. PRESENT(forces)) THEN
514 : ! print eigenvalues from fock-layer
515 8176 : CALL group%sum(evals)
516 8176 : IF (pao%iw_fockev > 0) THEN
517 2000 : DO iatom = 1, natoms
518 2000 : WRITE (pao%iw_fockev, *) "PAO| atom:", iatom, " fock evals around gap:", evals(:, iatom)
519 : END DO
520 500 : CALL m_flush(pao%iw_fockev)
521 : END IF
522 : ! print homo-lumo gap encountered by fock-layer
523 8176 : CALL group%min(gaps)
524 8176 : IF (pao%iw_gap > 0) THEN
525 2000 : DO iatom = 1, natoms
526 2000 : WRITE (pao%iw_gap, *) "PAO| atom:", iatom, " fock gap:", gaps(iatom)
527 : END DO
528 : END IF
529 : ! one-line summaries
530 8176 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| linpot regularization energy:", regu_energy
531 37804 : IF (pao%iw > 0) WRITE (pao%iw, "(A,E20.10,A,T71,I10)") " PAO| min_gap:", MINVAL(gaps), " for atom:", MINLOC(gaps)
532 : END IF
533 :
534 8210 : DEALLOCATE (gaps, evals)
535 8210 : CALL timestop(handle)
536 :
537 16420 : END SUBROUTINE pao_calc_U_linpot
538 :
539 : ! **************************************************************************************************
540 : !> \brief Internal routine, calculates terms in potential parametrization
541 : !> \param pao ...
542 : !> \param qs_env ...
543 : !> \param iatom ...
544 : !> \param V_blocks ...
545 : ! **************************************************************************************************
546 234 : SUBROUTINE linpot_calc_terms(pao, qs_env, iatom, V_blocks)
547 : TYPE(pao_env_type), POINTER :: pao
548 : TYPE(qs_environment_type), POINTER :: qs_env
549 : INTEGER, INTENT(IN) :: iatom
550 : REAL(dp), DIMENSION(:, :, :), INTENT(OUT) :: V_blocks
551 :
552 273 : SELECT CASE (pao%parameterization)
553 : CASE (pao_fock_param)
554 39 : CALL linpot_full_calc_terms(V_blocks)
555 : CASE (pao_rotinv_param)
556 195 : CALL linpot_rotinv_calc_terms(qs_env, iatom, V_blocks)
557 : CASE DEFAULT
558 234 : CPABORT("unknown parameterization")
559 : END SELECT
560 :
561 234 : END SUBROUTINE linpot_calc_terms
562 :
563 : ! **************************************************************************************************
564 : !> \brief Internal routine, calculates force contributions from potential parametrization
565 : !> \param pao ...
566 : !> \param qs_env ...
567 : !> \param iatom ...
568 : !> \param M_blocks ...
569 : !> \param forces ...
570 : ! **************************************************************************************************
571 34 : SUBROUTINE linpot_calc_forces(pao, qs_env, iatom, M_blocks, forces)
572 : TYPE(pao_env_type), POINTER :: pao
573 : TYPE(qs_environment_type), POINTER :: qs_env
574 : INTEGER, INTENT(IN) :: iatom
575 : REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: M_blocks
576 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: forces
577 :
578 66 : SELECT CASE (pao%parameterization)
579 : CASE (pao_fock_param)
580 : ! no force contributions
581 : CASE (pao_rotinv_param)
582 32 : CALL linpot_rotinv_calc_forces(qs_env, iatom, M_blocks, forces)
583 : CASE DEFAULT
584 34 : CPABORT("unknown parameterization")
585 : END SELECT
586 :
587 34 : END SUBROUTINE linpot_calc_forces
588 :
589 : ! **************************************************************************************************
590 : !> \brief Calculate initial guess for matrix_X
591 : !> \param pao ...
592 : !> \param qs_env ...
593 : ! **************************************************************************************************
594 34 : SUBROUTINE pao_param_initguess_linpot(pao, qs_env)
595 : TYPE(pao_env_type), POINTER :: pao
596 : TYPE(qs_environment_type), POINTER :: qs_env
597 :
598 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_initguess_linpot'
599 :
600 : INTEGER :: acol, arow, handle, i, iatom, j, k, n, &
601 : nterms
602 34 : INTEGER, DIMENSION(:), POINTER :: pri_basis_size
603 : LOGICAL :: found
604 : REAL(dp) :: w
605 34 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
606 34 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: S, S_evecs, S_inv
607 34 : REAL(dp), DIMENSION(:), POINTER :: V_guess_vec
608 34 : REAL(dp), DIMENSION(:, :), POINTER :: block_X, V_guess, V_terms
609 : TYPE(dbcsr_iterator_type) :: iter
610 :
611 34 : CALL timeset(routineN, handle)
612 :
613 34 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis_size)
614 :
615 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env,pri_basis_size) &
616 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)
617 : CALL dbcsr_iterator_start(iter, pao%matrix_X)
618 : DO WHILE (dbcsr_iterator_blocks_left(iter))
619 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
620 : iatom = arow; CPASSERT(arow == acol)
621 : CALL dbcsr_get_block_p(matrix=pao%matrix_V_terms, row=iatom, col=iatom, block=V_terms, found=found)
622 : CPASSERT(ASSOCIATED(V_terms))
623 : nterms = SIZE(V_terms, 2)
624 : IF (nterms == 0) CYCLE ! protect against corner-case of zero pao parameters
625 :
626 : ! guess initial potential
627 : N = pri_basis_size(iatom)
628 : ALLOCATE (V_guess_vec(n*n))
629 : V_guess(1:n, 1:n) => V_guess_vec
630 : CALL pao_guess_initial_potential(qs_env, iatom, V_guess)
631 :
632 : ! build overlap matrix
633 : ALLOCATE (S(nterms, nterms))
634 : S(:, :) = MATMUL(TRANSPOSE(V_terms), V_terms)
635 :
636 : ! diagonalize S
637 : ALLOCATE (S_evals(nterms), S_evecs(nterms, nterms))
638 : S_evecs(:, :) = S
639 : CALL diamat_all(S_evecs, S_evals)
640 :
641 : ! calculate Tikhonov regularized inverse
642 : ALLOCATE (S_inv(nterms, nterms))
643 : S_inv(:, :) = 0.0_dp
644 : DO k = 1, nterms
645 : w = S_evals(k)/(S_evals(k)**2 + pao%linpot_init_delta)
646 : DO i = 1, nterms
647 : DO j = 1, nterms
648 : S_inv(i, j) = S_inv(i, j) + w*S_evecs(i, k)*S_evecs(j, k)
649 : END DO
650 : END DO
651 : END DO
652 :
653 : ! perform fit
654 : block_X(:, 1) = MATMUL(MATMUL(S_inv, TRANSPOSE(V_terms)), V_guess_vec)
655 :
656 : ! clean up
657 : DEALLOCATE (V_guess_vec, S, S_evecs, S_evals, S_inv)
658 : END DO
659 : CALL dbcsr_iterator_stop(iter)
660 : !$OMP END PARALLEL
661 :
662 34 : CALL timestop(handle)
663 68 : END SUBROUTINE pao_param_initguess_linpot
664 :
665 24248 : END MODULE pao_param_linpot
|
