1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate MAO's and analyze wavefunctions 8!> \par History 9!> 03.2016 created [JGH] 10!> 12.2016 split into four modules [JGH] 11!> \author JGH 12! ************************************************************************************************** 13MODULE mao_methods 14 USE atomic_kind_types, ONLY: get_atomic_kind 15 USE basis_set_container_types, ONLY: add_basis_set_to_container 16 USE basis_set_types, ONLY: create_primitive_basis_set,& 17 get_gto_basis_set,& 18 gto_basis_set_p_type,& 19 gto_basis_set_type,& 20 write_gto_basis_set 21 USE cp_control_types, ONLY: dft_control_type 22 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 23 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 24 cp_dbcsr_plus_fm_fm_t,& 25 dbcsr_allocate_matrix_set 26 USE cp_fm_diag, ONLY: cp_fm_geeig 27 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 28 cp_fm_struct_release,& 29 cp_fm_struct_type 30 USE cp_fm_types, ONLY: cp_fm_create,& 31 cp_fm_release,& 32 cp_fm_type 33 USE cp_para_types, ONLY: cp_para_env_type 34 USE dbcsr_api, ONLY: & 35 dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_dot, dbcsr_get_block_p, & 36 dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, & 37 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, & 38 dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_set, dbcsr_type, & 39 dbcsr_type_no_symmetry 40 USE input_constants, ONLY: mao_basis_ext,& 41 mao_basis_orb,& 42 mao_basis_prim 43 USE iterate_matrix, ONLY: invert_Hotelling 44 USE kinds, ONLY: dp 45 USE kpoint_methods, ONLY: rskp_transform 46 USE kpoint_types, ONLY: get_kpoint_info,& 47 kpoint_type 48 USE lapack, ONLY: lapack_ssyev,& 49 lapack_ssygv 50 USE message_passing, ONLY: mp_sum 51 USE particle_types, ONLY: particle_type 52 USE qs_environment_types, ONLY: get_qs_env,& 53 qs_environment_type 54 USE qs_interactions, ONLY: init_interaction_radii_orb_basis 55 USE qs_kind_types, ONLY: get_qs_kind,& 56 qs_kind_type 57 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 58#include "./base/base_uses.f90" 59 60 IMPLICIT NONE 61 PRIVATE 62 63 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_methods' 64 65 TYPE mblocks 66 INTEGER :: n, ma 67 REAL(KIND=dp), DIMENSION(:, :), POINTER :: mat 68 REAL(KIND=dp), DIMENSION(:), POINTER :: eig 69 END TYPE mblocks 70 71 PUBLIC :: mao_initialization, mao_function, mao_function_gradient, mao_orthogonalization, & 72 mao_project_gradient, mao_scalar_product, mao_build_q, mao_basis_analysis, & 73 mao_reference_basis, calculate_p_gamma 74 75! ************************************************************************************************** 76 77CONTAINS 78 79! ************************************************************************************************** 80!> \brief ... 81!> \param mao_coef ... 82!> \param pmat ... 83!> \param smat ... 84!> \param eps1 ... 85! ************************************************************************************************** 86 SUBROUTINE mao_initialization(mao_coef, pmat, smat, eps1) 87 TYPE(dbcsr_type) :: mao_coef, pmat, smat 88 REAL(KIND=dp), INTENT(IN) :: eps1 89 90 CHARACTER(len=*), PARAMETER :: routineN = 'mao_initialization', & 91 routineP = moduleN//':'//routineN 92 93 INTEGER :: group, i, iatom, info, jatom, lwork, m, & 94 n, nblk 95 INTEGER, DIMENSION(:), POINTER :: mao_blk, row_blk, row_blk_sizes 96 LOGICAL :: found 97 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work 98 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat 99 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, pblock, sblock 100 TYPE(dbcsr_distribution_type) :: dbcsr_dist 101 TYPE(dbcsr_iterator_type) :: dbcsr_iter 102 TYPE(mblocks), ALLOCATABLE, DIMENSION(:) :: mbl 103 104 CALL dbcsr_get_info(mao_coef, nblkrows_total=nblk) 105 ALLOCATE (mbl(nblk)) 106 DO i = 1, nblk 107 NULLIFY (mbl(i)%mat, mbl(i)%eig) 108 END DO 109 110 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef) 111 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 112 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock) 113 CPASSERT(iatom == jatom) 114 m = SIZE(cblock, 2) 115 NULLIFY (pblock, sblock) 116 CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pblock, found=found) 117 CPASSERT(found) 118 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found) 119 CPASSERT(found) 120 n = SIZE(sblock, 1) 121 lwork = MAX(n*n, 100) 122 ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork)) 123 amat(1:n, 1:n) = pblock(1:n, 1:n) 124 bmat(1:n, 1:n) = sblock(1:n, 1:n) 125 info = 0 126 CALL lapack_ssygv(1, "V", "U", n, amat, n, bmat, n, w, work, lwork, info) 127 CPASSERT(info == 0) 128 ALLOCATE (mbl(iatom)%mat(n, n), mbl(iatom)%eig(n)) 129 mbl(iatom)%n = n 130 mbl(iatom)%ma = m 131 DO i = 1, n 132 mbl(iatom)%eig(i) = w(n - i + 1) 133 mbl(iatom)%mat(1:n, i) = amat(1:n, n - i + 1) 134 END DO 135 cblock(1:n, 1:m) = amat(1:n, n:n - m + 1:-1) 136 DEALLOCATE (amat, bmat, w, work) 137 END DO 138 CALL dbcsr_iterator_stop(dbcsr_iter) 139 140 IF (eps1 < 10.0_dp) THEN 141 CALL dbcsr_get_info(mao_coef, row_blk_size=row_blk_sizes, group=group) 142 ALLOCATE (row_blk(nblk), mao_blk(nblk)) 143 mao_blk = 0 144 row_blk = row_blk_sizes 145 DO iatom = 1, nblk 146 IF (ASSOCIATED(mbl(iatom)%mat)) THEN 147 n = mbl(iatom)%n 148 m = 0 149 DO i = 1, n 150 IF (mbl(iatom)%eig(i) < eps1) EXIT 151 m = i 152 END DO 153 m = MAX(m, mbl(iatom)%ma) 154 mbl(iatom)%ma = m 155 mao_blk(iatom) = m 156 END IF 157 END DO 158 CALL mp_sum(mao_blk, group) 159 CALL dbcsr_get_info(mao_coef, distribution=dbcsr_dist) 160 CALL dbcsr_release(mao_coef) 161 CALL dbcsr_create(mao_coef, name="MAO_COEF", dist=dbcsr_dist, & 162 matrix_type=dbcsr_type_no_symmetry, row_blk_size=row_blk, & 163 col_blk_size=mao_blk, nze=0) 164 CALL dbcsr_reserve_diag_blocks(matrix=mao_coef) 165 DEALLOCATE (mao_blk, row_blk) 166 ! 167 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef) 168 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 169 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock) 170 CPASSERT(iatom == jatom) 171 n = SIZE(cblock, 1) 172 m = SIZE(cblock, 2) 173 CPASSERT(n == mbl(iatom)%n .AND. m == mbl(iatom)%ma) 174 cblock(1:n, 1:m) = mbl(iatom)%mat(1:n, 1:m) 175 END DO 176 CALL dbcsr_iterator_stop(dbcsr_iter) 177 ! 178 END IF 179 180 CALL mao_orthogonalization(mao_coef, smat) 181 182 DO i = 1, nblk 183 IF (ASSOCIATED(mbl(i)%mat)) THEN 184 DEALLOCATE (mbl(i)%mat) 185 END IF 186 IF (ASSOCIATED(mbl(i)%eig)) THEN 187 DEALLOCATE (mbl(i)%eig) 188 END IF 189 END DO 190 DEALLOCATE (mbl) 191 192 END SUBROUTINE mao_initialization 193 194! ************************************************************************************************** 195!> \brief ... 196!> \param mao_coef ... 197!> \param fval ... 198!> \param qmat ... 199!> \param smat ... 200!> \param binv ... 201!> \param reuse ... 202! ************************************************************************************************** 203 SUBROUTINE mao_function(mao_coef, fval, qmat, smat, binv, reuse) 204 TYPE(dbcsr_type) :: mao_coef 205 REAL(KIND=dp), INTENT(OUT) :: fval 206 TYPE(dbcsr_type) :: qmat, smat, binv 207 LOGICAL, INTENT(IN) :: reuse 208 209 CHARACTER(len=*), PARAMETER :: routineN = 'mao_function', routineP = moduleN//':'//routineN 210 211 REAL(KIND=dp) :: convergence, threshold 212 TYPE(dbcsr_type) :: bmat, scmat, tmat 213 214 threshold = 1.e-8_dp 215 convergence = 1.e-6_dp 216 ! temp matrices 217 CALL dbcsr_create(scmat, template=mao_coef) 218 CALL dbcsr_create(bmat, template=binv) 219 CALL dbcsr_create(tmat, template=qmat) 220 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap 221 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat) 222 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat) 223 ! calculate inverse of B 224 CALL invert_Hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, & 225 norm_convergence=convergence, silent=.TRUE.) 226 ! calculate Binv*C and T=C(T)*Binv*C 227 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat) 228 CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat) 229 ! function = Tr(Q*T) 230 CALL dbcsr_dot(qmat, tmat, fval) 231 ! free temp matrices 232 CALL dbcsr_release(scmat) 233 CALL dbcsr_release(bmat) 234 CALL dbcsr_release(tmat) 235 236 END SUBROUTINE mao_function 237 238! ************************************************************************************************** 239!> \brief ... 240!> \param mao_coef ... 241!> \param fval ... 242!> \param mao_grad ... 243!> \param qmat ... 244!> \param smat ... 245!> \param binv ... 246!> \param reuse ... 247! ************************************************************************************************** 248 SUBROUTINE mao_function_gradient(mao_coef, fval, mao_grad, qmat, smat, binv, reuse) 249 TYPE(dbcsr_type) :: mao_coef 250 REAL(KIND=dp), INTENT(OUT) :: fval 251 TYPE(dbcsr_type) :: mao_grad, qmat, smat, binv 252 LOGICAL, INTENT(IN) :: reuse 253 254 CHARACTER(len=*), PARAMETER :: routineN = 'mao_function_gradient', & 255 routineP = moduleN//':'//routineN 256 257 REAL(KIND=dp) :: convergence, threshold 258 TYPE(dbcsr_type) :: bmat, scmat, t2mat, tmat 259 260 threshold = 1.e-8_dp 261 convergence = 1.e-6_dp 262 ! temp matrices 263 CALL dbcsr_create(scmat, template=mao_coef) 264 CALL dbcsr_create(bmat, template=binv) 265 CALL dbcsr_create(tmat, template=qmat) 266 CALL dbcsr_create(t2mat, template=scmat) 267 ! calculate B=C(T)*S*C matrix, S=(MAO,MAO) overlap 268 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef, 0.0_dp, scmat) 269 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef, scmat, 0.0_dp, bmat) 270 ! calculate inverse of B 271 CALL invert_Hotelling(binv, bmat, threshold, use_inv_as_guess=reuse, & 272 norm_convergence=convergence, silent=.TRUE.) 273 ! calculate R=C*Binv and T=C*Binv*C(T)=R*C(T) 274 CALL dbcsr_multiply("N", "N", 1.0_dp, mao_coef, binv, 0.0_dp, scmat) 275 CALL dbcsr_multiply("N", "T", 1.0_dp, scmat, mao_coef, 0.0_dp, tmat) 276 ! function = Tr(Q*T) 277 CALL dbcsr_dot(qmat, tmat, fval) 278 ! Gradient part 1: g = 2*Q*C*Binv = 2*Q*R 279 CALL dbcsr_multiply("N", "N", 2.0_dp, qmat, scmat, 0.0_dp, mao_grad, & 280 retain_sparsity=.TRUE.) 281 ! Gradient part 2: g = -2*S*T*X; X = Q*R 282 CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, scmat, 0.0_dp, t2mat) 283 CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, t2mat, 0.0_dp, scmat) 284 CALL dbcsr_multiply("N", "N", -2.0_dp, smat, scmat, 1.0_dp, mao_grad, & 285 retain_sparsity=.TRUE.) 286 ! free temp matrices 287 CALL dbcsr_release(scmat) 288 CALL dbcsr_release(bmat) 289 CALL dbcsr_release(tmat) 290 CALL dbcsr_release(t2mat) 291 292 CALL mao_project_gradient(mao_coef, mao_grad, smat) 293 294 END SUBROUTINE mao_function_gradient 295 296! ************************************************************************************************** 297!> \brief ... 298!> \param mao_coef ... 299!> \param smat ... 300! ************************************************************************************************** 301 SUBROUTINE mao_orthogonalization(mao_coef, smat) 302 TYPE(dbcsr_type) :: mao_coef, smat 303 304 CHARACTER(len=*), PARAMETER :: routineN = 'mao_orthogonalization', & 305 routineP = moduleN//':'//routineN 306 307 INTEGER :: i, iatom, info, jatom, lwork, m, n 308 LOGICAL :: found 309 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work 310 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat 311 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, sblock 312 TYPE(dbcsr_iterator_type) :: dbcsr_iter 313 314 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef) 315 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 316 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock) 317 CPASSERT(iatom == jatom) 318 m = SIZE(cblock, 2) 319 n = SIZE(cblock, 1) 320 NULLIFY (sblock) 321 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found) 322 CPASSERT(found) 323 lwork = MAX(n*n, 100) 324 ALLOCATE (amat(n, m), bmat(m, m), w(m), work(lwork)) 325 amat(1:n, 1:m) = MATMUL(sblock(1:n, 1:n), cblock(1:n, 1:m)) 326 bmat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), amat(1:n, 1:m)) 327 info = 0 328 CALL lapack_ssyev("V", "U", m, bmat, m, w, work, lwork, info) 329 CPASSERT(info == 0) 330 CPASSERT(ALL(w > 0.0_dp)) 331 w = 1.0_dp/SQRT(w) 332 DO i = 1, m 333 amat(1:m, i) = bmat(1:m, i)*w(i) 334 END DO 335 bmat(1:m, 1:m) = MATMUL(amat(1:m, 1:m), TRANSPOSE(bmat(1:m, 1:m))) 336 cblock(1:n, 1:m) = MATMUL(cblock(1:n, 1:m), bmat(1:m, 1:m)) 337 DEALLOCATE (amat, bmat, w, work) 338 END DO 339 CALL dbcsr_iterator_stop(dbcsr_iter) 340 341 END SUBROUTINE mao_orthogonalization 342 343! ************************************************************************************************** 344!> \brief ... 345!> \param mao_coef ... 346!> \param mao_grad ... 347!> \param smat ... 348! ************************************************************************************************** 349 SUBROUTINE mao_project_gradient(mao_coef, mao_grad, smat) 350 TYPE(dbcsr_type) :: mao_coef, mao_grad, smat 351 352 CHARACTER(len=*), PARAMETER :: routineN = 'mao_project_gradient', & 353 routineP = moduleN//':'//routineN 354 355 INTEGER :: iatom, jatom, m, n 356 LOGICAL :: found 357 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat 358 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cblock, gblock, sblock 359 TYPE(dbcsr_iterator_type) :: dbcsr_iter 360 361 CALL dbcsr_iterator_start(dbcsr_iter, mao_coef) 362 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 363 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, cblock) 364 CPASSERT(iatom == jatom) 365 m = SIZE(cblock, 2) 366 n = SIZE(cblock, 1) 367 NULLIFY (sblock) 368 CALL dbcsr_get_block_p(matrix=smat, row=iatom, col=jatom, block=sblock, found=found) 369 CPASSERT(found) 370 NULLIFY (gblock) 371 CALL dbcsr_get_block_p(matrix=mao_grad, row=iatom, col=jatom, block=gblock, found=found) 372 CPASSERT(found) 373 ALLOCATE (amat(m, m)) 374 amat(1:m, 1:m) = MATMUL(TRANSPOSE(cblock(1:n, 1:m)), MATMUL(sblock(1:n, 1:n), gblock(1:n, 1:m))) 375 gblock(1:n, 1:m) = gblock(1:n, 1:m) - MATMUL(cblock(1:n, 1:m), amat(1:m, 1:m)) 376 DEALLOCATE (amat) 377 END DO 378 CALL dbcsr_iterator_stop(dbcsr_iter) 379 380 END SUBROUTINE mao_project_gradient 381 382! ************************************************************************************************** 383!> \brief ... 384!> \param fmat1 ... 385!> \param fmat2 ... 386!> \return ... 387! ************************************************************************************************** 388 FUNCTION mao_scalar_product(fmat1, fmat2) RESULT(spro) 389 TYPE(dbcsr_type) :: fmat1, fmat2 390 REAL(KIND=dp) :: spro 391 392 CHARACTER(len=*), PARAMETER :: routineN = 'mao_scalar_product', & 393 routineP = moduleN//':'//routineN 394 395 INTEGER :: group, iatom, jatom, m, n 396 LOGICAL :: found 397 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ablock, bblock 398 TYPE(dbcsr_iterator_type) :: dbcsr_iter 399 400 spro = 0.0_dp 401 402 CALL dbcsr_iterator_start(dbcsr_iter, fmat1) 403 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 404 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, ablock) 405 CPASSERT(iatom == jatom) 406 m = SIZE(ablock, 2) 407 n = SIZE(ablock, 1) 408 CALL dbcsr_get_block_p(matrix=fmat2, row=iatom, col=jatom, block=bblock, found=found) 409 CPASSERT(found) 410 spro = spro + SUM(ablock(1:n, 1:m)*bblock(1:n, 1:m)) 411 END DO 412 CALL dbcsr_iterator_stop(dbcsr_iter) 413 414 CALL dbcsr_get_info(fmat1, group=group) 415 CALL mp_sum(spro, group) 416 417 END FUNCTION mao_scalar_product 418 419! ************************************************************************************************** 420!> \brief Calculate the density matrix at the Gamma point 421!> \param pmat ... 422!> \param ksmat ... 423!> \param smat ... 424!> \param kpoints Kpoint environment 425!> \param nmos Number of occupied orbitals 426!> \param occ Maximum occupation per orbital 427!> \par History 428!> 04.2016 created [JGH] 429! ************************************************************************************************** 430 SUBROUTINE calculate_p_gamma(pmat, ksmat, smat, kpoints, nmos, occ) 431 432 TYPE(dbcsr_type) :: pmat, ksmat, smat 433 TYPE(kpoint_type), POINTER :: kpoints 434 INTEGER, INTENT(IN) :: nmos 435 REAL(KIND=dp), INTENT(IN) :: occ 436 437 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_p_gamma', & 438 routineP = moduleN//':'//routineN 439 440 INTEGER :: norb 441 REAL(KIND=dp) :: de 442 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues 443 TYPE(cp_fm_struct_type), POINTER :: matrix_struct 444 TYPE(cp_fm_type), POINTER :: fmksmat, fmsmat, fmvec, fmwork 445 TYPE(dbcsr_type) :: tempmat 446 447 ! FM matrices 448 449 CALL dbcsr_get_info(smat, nfullrows_total=norb) 450 CALL cp_fm_struct_create(fmstruct=matrix_struct, context=kpoints%blacs_env_all, & 451 nrow_global=norb, ncol_global=norb) 452 CALL cp_fm_create(fmksmat, matrix_struct) 453 CALL cp_fm_create(fmsmat, matrix_struct) 454 CALL cp_fm_create(fmvec, matrix_struct) 455 CALL cp_fm_create(fmwork, matrix_struct) 456 ALLOCATE (eigenvalues(norb)) 457 458 ! DBCSR matrix 459 CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry) 460 461 ! transfer to FM 462 CALL dbcsr_desymmetrize(smat, tempmat) 463 CALL copy_dbcsr_to_fm(tempmat, fmsmat) 464 CALL dbcsr_desymmetrize(ksmat, tempmat) 465 CALL copy_dbcsr_to_fm(tempmat, fmksmat) 466 467 ! diagonalize 468 CALL cp_fm_geeig(fmksmat, fmsmat, fmvec, eigenvalues, fmwork) 469 de = eigenvalues(nmos + 1) - eigenvalues(nmos) 470 IF (de < 0.001_dp) THEN 471 CALL cp_warn(__LOCATION__, "MAO: No band gap at "// & 472 "Gamma point. MAO analysis not reliable.") 473 END IF 474 ! density matrix 475 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pmat, matrix_v=fmvec, ncol=nmos, alpha=occ) 476 477 DEALLOCATE (eigenvalues) 478 CALL dbcsr_release(tempmat) 479 CALL cp_fm_release(fmksmat) 480 CALL cp_fm_release(fmsmat) 481 CALL cp_fm_release(fmvec) 482 CALL cp_fm_release(fmwork) 483 CALL cp_fm_struct_release(matrix_struct) 484 485 END SUBROUTINE calculate_p_gamma 486 487! ************************************************************************************************** 488!> \brief Define the MAO reference basis set 489!> \param qs_env ... 490!> \param mao_basis ... 491!> \param mao_basis_set_list ... 492!> \param orb_basis_set_list ... 493!> \param iunit ... 494!> \param print_basis ... 495!> \par History 496!> 07.2016 created [JGH] 497! ************************************************************************************************** 498 SUBROUTINE mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, & 499 iunit, print_basis) 500 501 TYPE(qs_environment_type), POINTER :: qs_env 502 INTEGER, INTENT(IN) :: mao_basis 503 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list 504 INTEGER, INTENT(IN), OPTIONAL :: iunit 505 LOGICAL, INTENT(IN), OPTIONAL :: print_basis 506 507 CHARACTER(len=*), PARAMETER :: routineN = 'mao_reference_basis', & 508 routineP = moduleN//':'//routineN 509 510 INTEGER :: ikind, nbas, nkind, unit_nr 511 REAL(KIND=dp) :: eps_pgf_orb 512 TYPE(dft_control_type), POINTER :: dft_control 513 TYPE(gto_basis_set_type), POINTER :: basis_set, pbasis 514 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 515 TYPE(qs_kind_type), POINTER :: qs_kind 516 517 ! Reference basis set 518 CPASSERT(.NOT. ASSOCIATED(mao_basis_set_list)) 519 CPASSERT(.NOT. ASSOCIATED(orb_basis_set_list)) 520 521 ! options 522 IF (PRESENT(iunit)) THEN 523 unit_nr = iunit 524 ELSE 525 unit_nr = -1 526 END IF 527 528 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set) 529 nkind = SIZE(qs_kind_set) 530 ALLOCATE (mao_basis_set_list(nkind), orb_basis_set_list(nkind)) 531 DO ikind = 1, nkind 532 NULLIFY (mao_basis_set_list(ikind)%gto_basis_set) 533 NULLIFY (orb_basis_set_list(ikind)%gto_basis_set) 534 END DO 535 ! 536 DO ikind = 1, nkind 537 qs_kind => qs_kind_set(ikind) 538 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB") 539 IF (ASSOCIATED(basis_set)) orb_basis_set_list(ikind)%gto_basis_set => basis_set 540 END DO 541 ! 542 SELECT CASE (mao_basis) 543 CASE (mao_basis_orb) 544 DO ikind = 1, nkind 545 qs_kind => qs_kind_set(ikind) 546 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB") 547 IF (ASSOCIATED(basis_set)) mao_basis_set_list(ikind)%gto_basis_set => basis_set 548 END DO 549 CASE (mao_basis_prim) 550 DO ikind = 1, nkind 551 qs_kind => qs_kind_set(ikind) 552 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="ORB") 553 NULLIFY (pbasis) 554 IF (ASSOCIATED(basis_set)) THEN 555 CALL create_primitive_basis_set(basis_set, pbasis) 556 CALL get_qs_env(qs_env, dft_control=dft_control) 557 eps_pgf_orb = dft_control%qs_control%eps_pgf_orb 558 CALL init_interaction_radii_orb_basis(pbasis, eps_pgf_orb) 559 pbasis%kind_radius = basis_set%kind_radius 560 mao_basis_set_list(ikind)%gto_basis_set => pbasis 561 CALL add_basis_set_to_container(qs_kind%basis_sets, pbasis, "MAO") 562 END IF 563 END DO 564 CASE (mao_basis_ext) 565 DO ikind = 1, nkind 566 qs_kind => qs_kind_set(ikind) 567 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type="MAO") 568 IF (ASSOCIATED(basis_set)) THEN 569 basis_set%kind_radius = orb_basis_set_list(ikind)%gto_basis_set%kind_radius 570 mao_basis_set_list(ikind)%gto_basis_set => basis_set 571 END IF 572 END DO 573 CASE DEFAULT 574 CPABORT("Unknown option for MAO basis") 575 END SELECT 576 IF (unit_nr > 0) THEN 577 DO ikind = 1, nkind 578 IF (.NOT. ASSOCIATED(mao_basis_set_list(ikind)%gto_basis_set)) THEN 579 WRITE (UNIT=unit_nr, FMT="(T2,A,I4)") & 580 "WARNING: No MAO basis set associated with Kind ", ikind 581 ELSE 582 nbas = mao_basis_set_list(ikind)%gto_basis_set%nsgf 583 WRITE (UNIT=unit_nr, FMT="(T2,A,I4,T56,A,I10)") & 584 "MAO basis set Kind ", ikind, " Number of BSF:", nbas 585 END IF 586 END DO 587 END IF 588 589 IF (PRESENT(print_basis)) THEN 590 IF (print_basis) THEN 591 DO ikind = 1, nkind 592 basis_set => mao_basis_set_list(ikind)%gto_basis_set 593 CALL write_gto_basis_set(basis_set, unit_nr, "MAO REFERENCE BASIS") 594 END DO 595 END IF 596 END IF 597 598 END SUBROUTINE mao_reference_basis 599 600! ************************************************************************************************** 601!> \brief Analyze the MAO basis, projection on angular functions 602!> \param mao_coef ... 603!> \param matrix_smm ... 604!> \param mao_basis_set_list ... 605!> \param particle_set ... 606!> \param qs_kind_set ... 607!> \param unit_nr ... 608!> \param para_env ... 609!> \par History 610!> 07.2016 created [JGH] 611! ************************************************************************************************** 612 SUBROUTINE mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, & 613 qs_kind_set, unit_nr, para_env) 614 615 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, matrix_smm 616 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list 617 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 618 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 619 INTEGER, INTENT(IN) :: unit_nr 620 TYPE(cp_para_env_type), POINTER :: para_env 621 622 CHARACTER(len=*), PARAMETER :: routineN = 'mao_basis_analysis', & 623 routineP = moduleN//':'//routineN 624 625 CHARACTER(len=2) :: element_symbol 626 INTEGER :: ia, iab, iatom, ikind, iset, ishell, & 627 ispin, l, lmax, lshell, m, ma, na, & 628 natom, nspin 629 LOGICAL :: found 630 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cmask, vec1, vec2 631 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weight 632 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block, cmao 633 TYPE(gto_basis_set_type), POINTER :: basis_set 634 635 ! Analyze the MAO basis 636 IF (unit_nr > 0) THEN 637 WRITE (unit_nr, "(/,A)") " Analyze angular momentum character of MAOs " 638 WRITE (unit_nr, "(T7,A,T15,A,T20,A,T40,A,T50,A,T60,A,T70,A,T80,A)") & 639 "ATOM", "Spin", "MAO", "S", "P", "D", "F", "G" 640 END IF 641 lmax = 4 ! analyze up to g-functions 642 natom = SIZE(particle_set) 643 nspin = SIZE(mao_coef) 644 DO iatom = 1, natom 645 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 646 element_symbol=element_symbol, kind_number=ikind) 647 basis_set => mao_basis_set_list(ikind)%gto_basis_set 648 CALL get_qs_kind(qs_kind_set(ikind), mao=na) 649 CALL get_gto_basis_set(basis_set, nsgf=ma) 650 ALLOCATE (cmask(ma), vec1(ma), vec2(ma), weight(0:lmax, na)) 651 weight = 0.0_dp 652 CALL dbcsr_get_block_p(matrix=matrix_smm(1)%matrix, row=iatom, col=iatom, & 653 block=block, found=found) 654 DO ispin = 1, nspin 655 CALL dbcsr_get_block_p(matrix=mao_coef(ispin)%matrix, row=iatom, col=iatom, & 656 block=cmao, found=found) 657 IF (found) THEN 658 DO l = 0, lmax 659 cmask = 0.0_dp 660 iab = 0 661 DO iset = 1, basis_set%nset 662 DO ishell = 1, basis_set%nshell(iset) 663 lshell = basis_set%l(ishell, iset) 664 DO m = -lshell, lshell 665 iab = iab + 1 666 IF (l == lshell) cmask(iab) = 1.0_dp 667 END DO 668 END DO 669 END DO 670 DO ia = 1, na 671 vec1(1:ma) = cmask*cmao(1:ma, ia) 672 vec2(1:ma) = MATMUL(block, vec1) 673 weight(l, ia) = SUM(vec1(1:ma)*vec2(1:ma)) 674 END DO 675 END DO 676 END IF 677 CALL mp_sum(weight, para_env%group) 678 IF (unit_nr > 0) THEN 679 DO ia = 1, na 680 IF (ispin == 1 .AND. ia == 1) THEN 681 WRITE (unit_nr, "(i6,T9,A2,T17,i2,T20,i3,T31,5F10.4)") & 682 iatom, element_symbol, ispin, ia, weight(0:lmax, ia) 683 ELSE 684 WRITE (unit_nr, "(T17,i2,T20,i3,T31,5F10.4)") ispin, ia, weight(0:lmax, ia) 685 END IF 686 END DO 687 END IF 688 END DO 689 DEALLOCATE (cmask, weight, vec1, vec2) 690 END DO 691 END SUBROUTINE mao_basis_analysis 692 693! ************************************************************************************************** 694!> \brief Calculte the Q=APA(T) matrix, A=(MAO,ORB) overlap 695!> \param matrix_q ... 696!> \param matrix_p ... 697!> \param matrix_s ... 698!> \param matrix_smm ... 699!> \param matrix_smo ... 700!> \param smm_list ... 701!> \param electra ... 702!> \param eps_filter ... 703!> \param nimages ... 704!> \param kpoints ... 705!> \param matrix_ks ... 706!> \param sab_orb ... 707!> \par History 708!> 08.2016 created [JGH] 709! ************************************************************************************************** 710 SUBROUTINE mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, & 711 electra, eps_filter, nimages, kpoints, matrix_ks, sab_orb) 712 713 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_q 714 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s 715 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_smm, matrix_smo 716 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 717 POINTER :: smm_list 718 REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: electra 719 REAL(KIND=dp), INTENT(IN) :: eps_filter 720 INTEGER, INTENT(IN), OPTIONAL :: nimages 721 TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints 722 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 723 POINTER :: matrix_ks 724 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 725 OPTIONAL, POINTER :: sab_orb 726 727 CHARACTER(len=*), PARAMETER :: routineN = 'mao_build_q', routineP = moduleN//':'//routineN 728 729 INTEGER :: im, ispin, nim, nocc, norb, nspin 730 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 731 REAL(KIND=dp) :: elex, xkp(3) 732 TYPE(dbcsr_type) :: ksmat, pmat, smat, tmat 733 734 nim = 1 735 IF (PRESENT(nimages)) nim = nimages 736 IF (nim > 1) THEN 737 CPASSERT(PRESENT(kpoints)) 738 CPASSERT(PRESENT(matrix_ks)) 739 CPASSERT(PRESENT(sab_orb)) 740 END IF 741 742 ! Reference 743 nspin = SIZE(matrix_p, 1) 744 DO ispin = 1, nspin 745 electra(ispin) = 0.0_dp 746 DO im = 1, nim 747 CALL dbcsr_dot(matrix_p(ispin, im)%matrix, matrix_s(1, im)%matrix, elex) 748 electra(ispin) = electra(ispin) + elex 749 END DO 750 END DO 751 752 ! Q matrix 753 NULLIFY (matrix_q) 754 CALL dbcsr_allocate_matrix_set(matrix_q, nspin) 755 DO ispin = 1, nspin 756 ALLOCATE (matrix_q(ispin)%matrix) 757 CALL dbcsr_create(matrix_q(ispin)%matrix, template=matrix_smm(1)%matrix) 758 CALL cp_dbcsr_alloc_block_from_nbl(matrix_q(ispin)%matrix, smm_list) 759 END DO 760 ! temp matrix 761 CALL dbcsr_create(tmat, template=matrix_smo(1)%matrix, matrix_type=dbcsr_type_no_symmetry) 762 ! Q=APA(T) 763 DO ispin = 1, nspin 764 IF (nim == 1) THEN 765 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, matrix_p(ispin, 1)%matrix, & 766 0.0_dp, tmat, filter_eps=eps_filter) 767 CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, & 768 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter) 769 ELSE 770 ! k-points 771 CALL dbcsr_create(pmat, template=matrix_s(1, 1)%matrix) 772 CALL dbcsr_create(smat, template=matrix_s(1, 1)%matrix) 773 CALL dbcsr_create(ksmat, template=matrix_s(1, 1)%matrix) 774 CALL cp_dbcsr_alloc_block_from_nbl(pmat, sab_orb) 775 CALL cp_dbcsr_alloc_block_from_nbl(smat, sab_orb) 776 CALL cp_dbcsr_alloc_block_from_nbl(ksmat, sab_orb) 777 NULLIFY (cell_to_index) 778 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 779 ! calculate density matrix at gamma point 780 xkp = 0.0_dp 781 ! transform KS and S matrices to the gamma point 782 CALL dbcsr_set(ksmat, 0.0_dp) 783 CALL rskp_transform(rmatrix=ksmat, rsmat=matrix_ks, ispin=ispin, & 784 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb) 785 CALL dbcsr_set(smat, 0.0_dp) 786 CALL rskp_transform(rmatrix=smat, rsmat=matrix_s, ispin=1, & 787 xkp=xkp, cell_to_index=cell_to_index, sab_nl=sab_orb) 788 norb = NINT(electra(ispin)) 789 nocc = MOD(2, nspin) + 1 790 CALL calculate_p_gamma(pmat, ksmat, smat, kpoints, norb, REAL(nocc, KIND=dp)) 791 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smo(1)%matrix, pmat, & 792 0.0_dp, tmat, filter_eps=eps_filter) 793 CALL dbcsr_multiply("N", "T", 1.0_dp, tmat, matrix_smo(1)%matrix, & 794 0.0_dp, matrix_q(ispin)%matrix, filter_eps=eps_filter) 795 CALL dbcsr_release(pmat) 796 CALL dbcsr_release(smat) 797 CALL dbcsr_release(ksmat) 798 END IF 799 END DO 800 ! free temp matrix 801 CALL dbcsr_release(tmat) 802 803 END SUBROUTINE mao_build_q 804 805END MODULE mao_methods 806