1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate localized minimal basis 8!> \par History 9!> 12.2016 created [JGH] 10!> \author JGH 11! ************************************************************************************************** 12MODULE minbas_methods 13 USE cp_blacs_env, ONLY: cp_blacs_env_type 14 USE cp_control_types, ONLY: dft_control_type 15 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 16 copy_fm_to_dbcsr,& 17 dbcsr_allocate_matrix_set,& 18 dbcsr_deallocate_matrix_set 19 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale 20 USE cp_fm_diag, ONLY: choose_eigv_solver,& 21 cp_fm_power 22 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 23 cp_fm_struct_release,& 24 cp_fm_struct_type 25 USE cp_fm_types, ONLY: cp_fm_create,& 26 cp_fm_get_diag,& 27 cp_fm_release,& 28 cp_fm_to_fm_submat,& 29 cp_fm_type 30 USE cp_gemm_interface, ONLY: cp_gemm 31 USE cp_para_types, ONLY: cp_para_env_type 32 USE dbcsr_api, ONLY: & 33 dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_iterator_blocks_left, & 34 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 35 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type, & 36 dbcsr_type_no_symmetry 37 USE kinds, ONLY: dp 38 USE lapack, ONLY: lapack_ssyev 39 USE mao_basis, ONLY: mao_generate_basis 40 USE particle_methods, ONLY: get_particle_set 41 USE particle_types, ONLY: particle_type 42 USE qs_environment_types, ONLY: get_qs_env,& 43 qs_environment_type 44 USE qs_kind_types, ONLY: qs_kind_type 45 USE qs_ks_types, ONLY: get_ks_env,& 46 qs_ks_env_type 47 USE qs_mo_types, ONLY: get_mo_set,& 48 mo_set_p_type 49#include "./base/base_uses.f90" 50 51 IMPLICIT NONE 52 PRIVATE 53 54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_methods' 55 56 PUBLIC :: minbas_calculation 57 58! ************************************************************************************************** 59 60CONTAINS 61 62! ************************************************************************************************** 63!> \brief ... 64!> \param qs_env ... 65!> \param mos ... 66!> \param quambo ... 67!> \param mao ... 68!> \param iounit ... 69!> \param full_ortho ... 70!> \param eps_filter ... 71! ************************************************************************************************** 72 SUBROUTINE minbas_calculation(qs_env, mos, quambo, mao, iounit, full_ortho, eps_filter) 73 TYPE(qs_environment_type), POINTER :: qs_env 74 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 75 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: quambo 76 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 77 POINTER :: mao 78 INTEGER, INTENT(IN), OPTIONAL :: iounit 79 LOGICAL, INTENT(IN), OPTIONAL :: full_ortho 80 REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter 81 82 CHARACTER(len=*), PARAMETER :: routineN = 'minbas_calculation', & 83 routineP = moduleN//':'//routineN 84 85 INTEGER :: handle, homo, i, iab, ispin, nao, natom, & 86 ndep, nmao, nmo, nmx, np, np1, nspin, & 87 nvirt, unit_nr 88 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes 89 LOGICAL :: do_minbas, my_full_ortho 90 REAL(KIND=dp) :: my_eps_filter 91 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dval, dvalo, dvalv, eigval 92 TYPE(cp_blacs_env_type), POINTER :: blacs_env 93 TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_b, fm_struct_c, & 94 fm_struct_d, fm_struct_e 95 TYPE(cp_fm_type), POINTER :: fm1, fm2, fm3, fm4, fm5, fm6, fm_mos, & 96 fma, fmb, fmwork 97 TYPE(cp_para_env_type), POINTER :: para_env 98 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 99 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef 100 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s 101 TYPE(dbcsr_type) :: smao, sortho 102 TYPE(dbcsr_type), POINTER :: smat 103 TYPE(dft_control_type), POINTER :: dft_control 104 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 105 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 106 TYPE(qs_ks_env_type), POINTER :: ks_env 107 108 CALL timeset(routineN, handle) 109 110 IF (PRESENT(iounit)) THEN 111 unit_nr = iounit 112 ELSE 113 unit_nr = -1 114 END IF 115 116 IF (PRESENT(full_ortho)) THEN 117 my_full_ortho = full_ortho 118 ELSE 119 my_full_ortho = .FALSE. 120 END IF 121 122 IF (PRESENT(eps_filter)) THEN 123 my_eps_filter = eps_filter 124 ELSE 125 my_eps_filter = 1.0e-10_dp 126 END IF 127 128 CALL get_qs_env(qs_env, dft_control=dft_control) 129 nspin = dft_control%nspins 130 131 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) 132 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom) 133 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist) 134 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) 135 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes) 136 CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes) 137 nmao = SUM(col_blk_sizes) 138 ! check if MAOs have been specified 139 DO iab = 1, natom 140 IF (col_blk_sizes(iab) < 0) & 141 CPABORT("Number of MAOs has to be specified in KIND section for all elements") 142 END DO 143 CALL get_mo_set(mo_set=mos(1)%mo_set, nao=nao, nmo=nmo) 144 145 IF (unit_nr > 0) THEN 146 WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Atomic Basis Set Functions :', nao 147 WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Minimal Basis Set Functions :', nmao 148 IF (nspin == 1) THEN 149 WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Molecular Orbitals available :', nmo 150 ELSE 151 DO ispin = 1, nspin 152 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nmo=nmx) 153 WRITE (unit_nr, '(T2,A,i2,T71,I10)') & 154 'Total Number of Molecular Orbitals available for Spin ', ispin, nmx 155 END DO 156 END IF 157 END IF 158 CPASSERT(nmao <= nao) 159 DO ispin = 1, nspin 160 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nmo=nmx) 161 IF (nmx /= nmo) EXIT 162 END DO 163 do_minbas = .TRUE. 164 IF (nmao > nmo) THEN 165 IF (unit_nr > 0) THEN 166 WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible' 167 END IF 168 do_minbas = .FALSE. 169 ELSEIF (nmo /= nmx) THEN 170 IF (unit_nr > 0) THEN 171 WRITE (unit_nr, '(T2,A)') 'Different Number of Alpha and Beta MOs' 172 WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible' 173 END IF 174 do_minbas = .FALSE. 175 ELSE 176 IF (nao > nmo) THEN 177 IF (unit_nr > 0) THEN 178 WRITE (unit_nr, '(T2,A)') 'WARNING: Only a subset of MOs is available: Analysis depends on MOs' 179 END IF 180 END IF 181 END IF 182 183 IF (do_minbas) THEN 184 ! initialize QUAMBOs 185 NULLIFY (quambo) 186 CALL dbcsr_allocate_matrix_set(quambo, nspin) 187 DO ispin = 1, nspin 188 ! coeficients 189 ALLOCATE (quambo(ispin)%matrix) 190 CALL dbcsr_create(matrix=quambo(ispin)%matrix, & 191 name="QUAMBO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, & 192 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 193 END DO 194 195 ! initialize MAOs 196 ! optimize MAOs (mao_coef is allocated in the routine) 197 CALL mao_generate_basis(qs_env, mao_coef) 198 199 ! sortho (nmao x nmao) 200 CALL dbcsr_create(sortho, name="SORTHO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, & 201 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 202 CALL dbcsr_reserve_diag_blocks(matrix=sortho) 203 204 DEALLOCATE (row_blk_sizes, col_blk_sizes) 205 206 ! temporary FM matrices 207 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) 208 NULLIFY (fm_struct_a, fm_struct_b) 209 CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, & 210 para_env=para_env, context=blacs_env) 211 CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmo, ncol_global=nmao, & 212 para_env=para_env, context=blacs_env) 213 CALL cp_fm_create(fm1, fm_struct_a) 214 CALL cp_fm_create(fm2, fm_struct_b) 215 CALL cp_fm_create(fma, fm_struct_b) 216 CALL cp_fm_create(fmb, fm_struct_b) 217 218 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s) 219 smat => matrix_s(1, 1)%matrix 220 DO ispin = 1, nspin 221 222 ! SMAO = Overlap*MAOs 223 CALL dbcsr_create(smao, name="S*MAO", template=mao_coef(1)%matrix) 224 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef(ispin)%matrix, 0.0_dp, smao) 225 ! a(nj)* = C(vn)(T) * SMAO(vj) 226 CALL copy_dbcsr_to_fm(smao, fm1) 227 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=fm_mos) 228 CALL cp_gemm("T", "N", nmo, nmao, nao, 1.0_dp, fm_mos, fm1, 0.0_dp, fm2) 229 CALL dbcsr_release(smao) 230 CALL get_mo_set(mo_set=mos(ispin)%mo_set, homo=homo) 231 IF (unit_nr > 0) THEN 232 WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Occupied Valence Set', 'Spin ', ispin, homo 233 END IF 234 nvirt = nmo - homo 235 NULLIFY (fm_struct_c) 236 CALL cp_fm_struct_create(fm_struct_c, nrow_global=nvirt, ncol_global=nvirt, & 237 para_env=para_env, context=blacs_env) 238 CALL cp_fm_create(fm3, fm_struct_c) 239 CALL cp_fm_create(fm4, fm_struct_c) 240 ! B(vw) = a(vj)* * a(wj)* 241 CALL cp_gemm("N", "T", nvirt, nvirt, nmao, 1.0_dp, fm2, fm2, 0.0_dp, fm3, & 242 a_first_row=homo + 1, b_first_row=homo + 1) 243 ALLOCATE (eigval(nvirt)) 244 CALL choose_eigv_solver(fm3, fm4, eigval) 245 ! SVD(B) -> select p largest eigenvalues and vectors 246 np = nmao - homo 247 np1 = nvirt - np + 1 248 IF (unit_nr > 0) THEN 249 WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Virtual Valence Set', 'Spin ', ispin, np 250 END IF 251 ! R(vw) = SUM_p T(vp)*T(wp) 252 CALL cp_gemm("N", "T", nvirt, nvirt, np, 1.0_dp, fm4, fm4, 0.0_dp, fm3, & 253 a_first_col=np1, b_first_col=np1) 254 ! 255 ALLOCATE (dval(nmao), dvalo(nmao), dvalv(nmao)) 256 NULLIFY (fm_struct_d) 257 CALL cp_fm_struct_create(fm_struct_d, nrow_global=nvirt, ncol_global=nmao, & 258 para_env=para_env, context=blacs_env) 259 CALL cp_fm_create(fm5, fm_struct_d) 260 NULLIFY (fm_struct_e) 261 CALL cp_fm_struct_create(fm_struct_e, nrow_global=nmao, ncol_global=nmao, & 262 para_env=para_env, context=blacs_env) 263 CALL cp_fm_create(fm6, fm_struct_e) 264 ! D(j) = SUM_n (a(nj)*)^2 + SUM_vw R(vw) * a(vj)* * a(wj)* 265 CALL cp_gemm("N", "N", nvirt, nmao, nvirt, 1.0_dp, fm3, fm2, 0.0_dp, fm5, & 266 b_first_row=homo + 1) 267 CALL cp_gemm("T", "N", nmao, nmao, nvirt, 1.0_dp, fm2, fm5, 0.0_dp, fm6, & 268 a_first_row=homo + 1) 269 CALL cp_fm_get_diag(fm6, dvalv(1:nmao)) 270 CALL cp_gemm("T", "N", nmao, nmao, homo, 1.0_dp, fm2, fm2, 0.0_dp, fm6) 271 CALL cp_fm_get_diag(fm6, dvalo(1:nmao)) 272 DO i = 1, nmao 273 dval(i) = 1.0_dp/SQRT(dvalo(i) + dvalv(i)) 274 END DO 275 ! scale intermediate expansion 276 CALL cp_fm_to_fm_submat(fm2, fma, homo, nmao, 1, 1, 1, 1) 277 CALL cp_fm_to_fm_submat(fm5, fma, nvirt, nmao, 1, 1, homo + 1, 1) 278 CALL cp_fm_column_scale(fma, dval) 279 ! Orthogonalization 280 CALL cp_fm_create(fmwork, fm_struct_e) 281 CALL cp_gemm("T", "N", nmao, nmao, nmo, 1.0_dp, fma, fma, 0.0_dp, fm6) 282 IF (my_full_ortho) THEN 283 ! full orthogonalization 284 CALL cp_fm_power(fm6, fmwork, -0.5_dp, 1.0e-12_dp, ndep) 285 IF (ndep > 0 .AND. unit_nr > 0) THEN 286 WRITE (unit_nr, '(T2,A,T71,I10)') 'Warning: linear dependent basis ', ndep 287 END IF 288 CALL cp_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb) 289 ELSE 290 ! orthogonalize on-atom blocks 291 CALL copy_fm_to_dbcsr(fm6, sortho, keep_sparsity=.TRUE.) 292 CALL diag_sqrt_invert(sortho) 293 CALL copy_dbcsr_to_fm(sortho, fm6) 294 CALL cp_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb) 295 END IF 296 ! store as QUAMBO 297 CALL cp_gemm("N", "N", nao, nmao, nmo, 1.0_dp, fm_mos, fmb, 0.0_dp, fm1) 298 CALL copy_fm_to_dbcsr(fm1, quambo(ispin)%matrix) 299 CALL dbcsr_filter(quambo(ispin)%matrix, my_eps_filter) 300 ! 301 DEALLOCATE (eigval, dval, dvalo, dvalv) 302 CALL cp_fm_release(fm3) 303 CALL cp_fm_release(fm4) 304 CALL cp_fm_release(fm5) 305 CALL cp_fm_release(fm6) 306 CALL cp_fm_release(fmwork) 307 CALL cp_fm_struct_release(fm_struct_c) 308 CALL cp_fm_struct_release(fm_struct_d) 309 CALL cp_fm_struct_release(fm_struct_e) 310 311 END DO 312 313 ! clean up 314 CALL cp_fm_release(fm1) 315 CALL cp_fm_release(fm2) 316 CALL cp_fm_release(fma) 317 CALL cp_fm_release(fmb) 318 CALL cp_fm_struct_release(fm_struct_a) 319 CALL cp_fm_struct_release(fm_struct_b) 320 CALL dbcsr_release(sortho) 321 322 ! return MAOs if requested 323 IF (PRESENT(mao)) THEN 324 mao => mao_coef 325 ELSE 326 CALL dbcsr_deallocate_matrix_set(mao_coef) 327 END IF 328 329 ELSE 330 NULLIFY (quambo) 331 END IF 332 333 CALL timestop(handle) 334 335 END SUBROUTINE minbas_calculation 336 337! ************************************************************************************************** 338!> \brief ... 339!> \param sortho ... 340! ************************************************************************************************** 341 SUBROUTINE diag_sqrt_invert(sortho) 342 TYPE(dbcsr_type) :: sortho 343 344 CHARACTER(len=*), PARAMETER :: routineN = 'diag_sqrt_invert', & 345 routineP = moduleN//':'//routineN 346 347 INTEGER :: i, iatom, info, jatom, lwork, n 348 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work 349 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat 350 REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock 351 TYPE(dbcsr_iterator_type) :: dbcsr_iter 352 353 CALL dbcsr_iterator_start(dbcsr_iter, sortho) 354 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 355 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, sblock) 356 CPASSERT(iatom == jatom) 357 n = SIZE(sblock, 1) 358 lwork = MAX(n*n, 100) 359 ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork)) 360 amat(1:n, 1:n) = sblock(1:n, 1:n) 361 info = 0 362 CALL lapack_ssyev("V", "U", n, amat, n, w, work, lwork, info) 363 CPASSERT(info == 0) 364 w(1:n) = 1._dp/SQRT(w(1:n)) 365 DO i = 1, n 366 bmat(1:n, i) = amat(1:n, i)*w(i) 367 END DO 368 sblock(1:n, 1:n) = MATMUL(amat, TRANSPOSE(bmat)) 369 DEALLOCATE (amat, bmat, w, work) 370 END DO 371 CALL dbcsr_iterator_stop(dbcsr_iter) 372 373 END SUBROUTINE diag_sqrt_invert 374 375END MODULE minbas_methods 376