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