1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of Overlap and Hamiltonian matrices in DFTB 8!> \author JGH 9! ************************************************************************************************** 10MODULE qs_dftb_matrices 11 USE atomic_kind_types, ONLY: atomic_kind_type,& 12 get_atomic_kind,& 13 get_atomic_kind_set 14 USE atprop_types, ONLY: atprop_array_init,& 15 atprop_type 16 USE block_p_types, ONLY: block_p_type 17 USE cp_control_types, ONLY: dft_control_type,& 18 dftb_control_type 19 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 20 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 21 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 22 USE cp_log_handling, ONLY: cp_get_default_logger,& 23 cp_logger_type 24 USE cp_output_handling, ONLY: cp_p_file,& 25 cp_print_key_finished_output,& 26 cp_print_key_should_output,& 27 cp_print_key_unit_nr 28 USE cp_para_types, ONLY: cp_para_env_type 29 USE dbcsr_api, ONLY: & 30 dbcsr_add, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, & 31 dbcsr_distribution_type, dbcsr_dot, dbcsr_finalize, dbcsr_get_block_p, dbcsr_multiply, & 32 dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_symmetric 33 USE efield_tb_methods, ONLY: efield_tb_matrix 34 USE input_section_types, ONLY: section_vals_get_subs_vals,& 35 section_vals_type,& 36 section_vals_val_get 37 USE kinds, ONLY: default_string_length,& 38 dp 39 USE kpoint_types, ONLY: get_kpoint_info,& 40 kpoint_type 41 USE message_passing, ONLY: mp_sum 42 USE mulliken, ONLY: mulliken_charges 43 USE particle_methods, ONLY: get_particle_set 44 USE particle_types, ONLY: particle_type 45 USE qs_dftb_coulomb, ONLY: build_dftb_coulomb 46 USE qs_dftb_types, ONLY: qs_dftb_atom_type,& 47 qs_dftb_pairpot_type 48 USE qs_dftb_utils, ONLY: compute_block_sk,& 49 get_dftb_atom_param,& 50 iptr,& 51 urep_egr 52 USE qs_energy_types, ONLY: qs_energy_type 53 USE qs_environment_types, ONLY: get_qs_env,& 54 qs_environment_type 55 USE qs_force_types, ONLY: qs_force_type 56 USE qs_kind_types, ONLY: get_qs_kind,& 57 get_qs_kind_set,& 58 qs_kind_type 59 USE qs_ks_types, ONLY: get_ks_env,& 60 qs_ks_env_type,& 61 set_ks_env 62 USE qs_mo_types, ONLY: get_mo_set,& 63 mo_set_p_type 64 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 65 neighbor_list_iterate,& 66 neighbor_list_iterator_create,& 67 neighbor_list_iterator_p_type,& 68 neighbor_list_iterator_release,& 69 neighbor_list_set_p_type 70 USE qs_rho_types, ONLY: qs_rho_get,& 71 qs_rho_type 72 USE virial_methods, ONLY: virial_pair_force 73 USE virial_types, ONLY: virial_type 74#include "./base/base_uses.f90" 75 76 IMPLICIT NONE 77 78 INTEGER, DIMENSION(16), PARAMETER :: orbptr = (/0, 1, 1, 1, & 79 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3/) 80 81 PRIVATE 82 83 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_matrices' 84 85 PUBLIC :: build_dftb_matrices, build_dftb_ks_matrix, build_dftb_overlap 86 87CONTAINS 88 89! ************************************************************************************************** 90!> \brief ... 91!> \param qs_env ... 92!> \param para_env ... 93!> \param calculate_forces ... 94! ************************************************************************************************** 95 SUBROUTINE build_dftb_matrices(qs_env, para_env, calculate_forces) 96 97 TYPE(qs_environment_type), POINTER :: qs_env 98 TYPE(cp_para_env_type), POINTER :: para_env 99 LOGICAL, INTENT(IN) :: calculate_forces 100 101 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dftb_matrices', & 102 routineP = moduleN//':'//routineN 103 104 INTEGER :: after, atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, irow, iw, jatom, & 105 jkind, l1, l2, la, lb, llm, lmaxi, lmaxj, m, n1, n2, n_urpoly, natom, natorb_a, natorb_b, & 106 nderivatives, ngrd, ngrdcut, nimg, nkind, spdim 107 INTEGER, DIMENSION(3) :: cell 108 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 109 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 110 LOGICAL :: defined, found, omit_headers, use_virial 111 REAL(KIND=dp) :: ddr, dgrd, dr, erep, erepij, f0, f1, & 112 foab, fow, s_cut, urep_cut 113 REAL(KIND=dp), DIMENSION(0:3) :: eta_a, eta_b, skself 114 REAL(KIND=dp), DIMENSION(10) :: urep 115 REAL(KIND=dp), DIMENSION(2) :: surr 116 REAL(KIND=dp), DIMENSION(3) :: drij, force_ab, force_rr, force_w, rij, & 117 srep 118 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dfblock, dsblock, fblock, fmatij, & 119 fmatji, pblock, sblock, scoeff, & 120 smatij, smatji, spxr, wblock 121 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 122 TYPE(atprop_type), POINTER :: atprop 123 TYPE(block_p_type), DIMENSION(2:4) :: dsblocks 124 TYPE(cp_logger_type), POINTER :: logger 125 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_w 126 TYPE(dft_control_type), POINTER :: dft_control 127 TYPE(dftb_control_type), POINTER :: dftb_control 128 TYPE(kpoint_type), POINTER :: kpoints 129 TYPE(neighbor_list_iterator_p_type), & 130 DIMENSION(:), POINTER :: nl_iterator 131 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 132 POINTER :: sab_orb 133 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 134 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b 135 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & 136 POINTER :: dftb_potential 137 TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji 138 TYPE(qs_energy_type), POINTER :: energy 139 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 140 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 141 TYPE(qs_ks_env_type), POINTER :: ks_env 142 TYPE(qs_rho_type), POINTER :: rho 143 TYPE(virial_type), POINTER :: virial 144 145 CALL timeset(routineN, handle) 146 147 ! set pointers 148 iptr = 0 149 DO la = 0, 3 150 DO lb = 0, 3 151 llm = 0 152 DO l1 = 0, MAX(la, lb) 153 DO l2 = 0, MIN(l1, la, lb) 154 DO m = 0, l2 155 llm = llm + 1 156 iptr(l1, l2, m, la, lb) = llm 157 END DO 158 END DO 159 END DO 160 END DO 161 END DO 162 163 NULLIFY (logger, virial, atprop) 164 logger => cp_get_default_logger() 165 166 NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, & 167 qs_kind_set, sab_orb, ks_env) 168 169 CALL get_qs_env(qs_env=qs_env, & 170 energy=energy, & 171 atomic_kind_set=atomic_kind_set, & 172 qs_kind_set=qs_kind_set, & 173 matrix_h_kp=matrix_h, & 174 matrix_s_kp=matrix_s, & 175 atprop=atprop, & 176 dft_control=dft_control, & 177 ks_env=ks_env) 178 179 dftb_control => dft_control%qs_control%dftb_control 180 nimg = dft_control%nimages 181 ! Allocate the overlap and Hamiltonian matrix 182 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb) 183 nderivatives = 0 184 IF (dftb_control%self_consistent .AND. calculate_forces) nderivatives = 1 185 CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s, "OVERLAP", sab_orb) 186 CALL setup_matrices2(qs_env, 0, nimg, matrix_h, "CORE HAMILTONIAN", sab_orb) 187 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s) 188 CALL set_ks_env(ks_env, matrix_h_kp=matrix_h) 189 190 NULLIFY (dftb_potential) 191 CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential) 192 NULLIFY (particle_set) 193 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set) 194 195 IF (calculate_forces) THEN 196 NULLIFY (rho, force, matrix_w) 197 CALL get_qs_env(qs_env=qs_env, & 198 rho=rho, & 199 matrix_w_kp=matrix_w, & 200 virial=virial, & 201 force=force) 202 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 203 204 IF (SIZE(matrix_p, 1) == 2) THEN 205 DO img = 1, nimg 206 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 207 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 208 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, & 209 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 210 END DO 211 END IF 212 natom = SIZE(particle_set) 213 ALLOCATE (atom_of_kind(natom)) 214 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 215 atom_of_kind=atom_of_kind) 216 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 217 END IF 218 ! atomic energy decomposition 219 IF (atprop%energy) THEN 220 natom = SIZE(particle_set) 221 CALL atprop_array_init(atprop%atecc, natom) 222 END IF 223 224 NULLIFY (cell_to_index) 225 IF (nimg > 1) THEN 226 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 227 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 228 END IF 229 230 erep = 0._dp 231 232 nkind = SIZE(atomic_kind_set) 233 234 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 235 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 236 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 237 iatom=iatom, jatom=jatom, r=rij, cell=cell) 238 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom) 239 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a) 240 CALL get_dftb_atom_param(dftb_kind_a, & 241 defined=defined, lmax=lmaxi, skself=skself, & 242 eta=eta_a, natorb=natorb_a) 243 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 244 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b) 245 CALL get_dftb_atom_param(dftb_kind_b, & 246 defined=defined, lmax=lmaxj, eta=eta_b, natorb=natorb_b) 247 248 IF (.NOT. defined .OR. natorb_b < 1) CYCLE 249 250 ! retrieve information on F and S matrix 251 dftb_param_ij => dftb_potential(ikind, jkind) 252 dftb_param_ji => dftb_potential(jkind, ikind) 253 ! assume table size and type is symmetric 254 ngrd = dftb_param_ij%ngrd 255 ngrdcut = dftb_param_ij%ngrdcut 256 dgrd = dftb_param_ij%dgrd 257 ddr = dgrd*0.1_dp 258 CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm) 259 llm = dftb_param_ij%llm 260 fmatij => dftb_param_ij%fmat 261 smatij => dftb_param_ij%smat 262 fmatji => dftb_param_ji%fmat 263 smatji => dftb_param_ji%smat 264 ! repulsive pair potential 265 n_urpoly = dftb_param_ij%n_urpoly 266 urep_cut = dftb_param_ij%urep_cut 267 urep = dftb_param_ij%urep 268 spxr => dftb_param_ij%spxr 269 scoeff => dftb_param_ij%scoeff 270 spdim = dftb_param_ij%spdim 271 s_cut = dftb_param_ij%s_cut 272 srep = dftb_param_ij%srep 273 surr = dftb_param_ij%surr 274 275 dr = SQRT(SUM(rij(:)**2)) 276 IF (NINT(dr/dgrd) <= ngrdcut) THEN 277 278 IF (nimg == 1) THEN 279 ic = 1 280 ELSE 281 ic = cell_to_index(cell(1), cell(2), cell(3)) 282 CPASSERT(ic > 0) 283 END IF 284 285 icol = MAX(iatom, jatom) 286 irow = MIN(iatom, jatom) 287 NULLIFY (sblock, fblock) 288 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, & 289 row=irow, col=icol, BLOCK=sblock, found=found) 290 CPASSERT(found) 291 CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, & 292 row=irow, col=icol, BLOCK=fblock, found=found) 293 CPASSERT(found) 294 295 IF (calculate_forces) THEN 296 NULLIFY (pblock) 297 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, & 298 row=irow, col=icol, block=pblock, found=found) 299 CPASSERT(ASSOCIATED(pblock)) 300 NULLIFY (wblock) 301 CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, & 302 row=irow, col=icol, block=wblock, found=found) 303 CPASSERT(ASSOCIATED(wblock)) 304 IF (dftb_control%self_consistent) THEN 305 DO i = 2, 4 306 NULLIFY (dsblocks(i)%block) 307 CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, & 308 row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found) 309 CPASSERT(found) 310 END DO 311 END IF 312 END IF 313 314 IF (iatom == jatom .AND. dr < 0.001_dp) THEN 315 ! diagonal block 316 DO i = 1, natorb_a 317 sblock(i, i) = sblock(i, i) + 1._dp 318 fblock(i, i) = fblock(i, i) + skself(orbptr(i)) 319 END DO 320 ELSE 321 ! off-diagonal block 322 CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, & 323 llm, lmaxi, lmaxj, irow, iatom) 324 CALL compute_block_sk(fblock, fmatij, fmatji, rij, ngrd, ngrdcut, dgrd, & 325 llm, lmaxi, lmaxj, irow, iatom) 326 IF (calculate_forces) THEN 327 force_ab = 0._dp 328 force_w = 0._dp 329 n1 = SIZE(fblock, 1) 330 n2 = SIZE(fblock, 2) 331 ! make sure that displacement is in the correct direction depending on the position 332 ! of the block (upper or lower triangle) 333 f0 = 1.0_dp 334 IF (irow == iatom) f0 = -1.0_dp 335 336 ALLOCATE (dfblock(n1, n2), dsblock(n1, n2)) 337 338 DO i = 1, 3 339 drij = rij 340 dfblock = 0._dp; dsblock = 0._dp 341 342 drij(i) = rij(i) - ddr*f0 343 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 344 llm, lmaxi, lmaxj, irow, iatom) 345 CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, & 346 llm, lmaxi, lmaxj, irow, iatom) 347 348 dsblock = -dsblock 349 dfblock = -dfblock 350 351 drij(i) = rij(i) + ddr*f0 352 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 353 llm, lmaxi, lmaxj, irow, iatom) 354 CALL compute_block_sk(dfblock, fmatij, fmatji, drij, ngrd, ngrdcut, dgrd, & 355 llm, lmaxi, lmaxj, irow, iatom) 356 357 dfblock = dfblock/(2.0_dp*ddr) 358 dsblock = dsblock/(2.0_dp*ddr) 359 360 foab = 2.0_dp*SUM(dfblock*pblock) 361 fow = -2.0_dp*SUM(dsblock*wblock) 362 363 force_ab(i) = force_ab(i) + foab 364 force_w(i) = force_w(i) + fow 365 IF (dftb_control%self_consistent) THEN 366 CPASSERT(ASSOCIATED(dsblocks(i + 1)%block)) 367 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock 368 END IF 369 ENDDO 370 IF (use_virial) THEN 371!deb iatom=jatom f0=0.5*f0 372 IF (iatom == jatom) f0 = 0.5_dp*f0 373!deb 374 CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij) 375 CALL virial_pair_force(virial%pv_virial, -f0, force_w, rij) 376 IF (atprop%stress) THEN 377 f1 = 0.5_dp*f0 378 CALL virial_pair_force(atprop%atstress(:, :, iatom), -f1, force_ab, rij) 379 CALL virial_pair_force(atprop%atstress(:, :, iatom), -f1, force_w, rij) 380 CALL virial_pair_force(atprop%atstress(:, :, jatom), -f1, force_ab, rij) 381 CALL virial_pair_force(atprop%atstress(:, :, jatom), -f1, force_w, rij) 382 END IF 383 END IF 384 DEALLOCATE (dfblock, dsblock) 385 END IF 386 END IF 387 388 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN 389 atom_a = atom_of_kind(iatom) 390 atom_b = atom_of_kind(jatom) 391 IF (irow == iatom) force_ab = -force_ab 392 IF (irow == iatom) force_w = -force_w 393 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:) 394 force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:) 395 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - force_w(:) 396 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + force_w(:) 397 END IF 398 399 END IF 400 401 ! repulsive potential 402 IF ((dr <= urep_cut .OR. spdim > 0) .AND. dr > 0.001_dp) THEN 403 erepij = 0._dp 404 CALL urep_egr(rij, dr, erepij, force_rr, & 405 n_urpoly, urep, spdim, s_cut, srep, spxr, scoeff, surr, calculate_forces) 406 erep = erep + erepij 407 IF (atprop%energy) THEN 408 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*erepij 409 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*erepij 410 END IF 411 IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN 412 atom_a = atom_of_kind(iatom) 413 atom_b = atom_of_kind(jatom) 414 force(ikind)%repulsive(:, atom_a) = & 415 force(ikind)%repulsive(:, atom_a) - force_rr(:) 416 force(jkind)%repulsive(:, atom_b) = & 417 force(jkind)%repulsive(:, atom_b) + force_rr(:) 418 IF (use_virial) THEN 419 f0 = -1.0_dp 420 IF (iatom == jatom) f0 = -0.5_dp 421 CALL virial_pair_force(virial%pv_virial, f0, force_rr, rij) 422 IF (atprop%stress) THEN 423 CALL virial_pair_force(atprop%atstress(:, :, iatom), f0*0.5_dp, force_rr, rij) 424 CALL virial_pair_force(atprop%atstress(:, :, jatom), f0*0.5_dp, force_rr, rij) 425 END IF 426 END IF 427 END IF 428 END IF 429 430 END DO 431 CALL neighbor_list_iterator_release(nl_iterator) 432 433 DO i = 1, SIZE(matrix_s, 1) 434 DO img = 1, nimg 435 CALL dbcsr_finalize(matrix_s(i, img)%matrix) 436 END DO 437 ENDDO 438 DO i = 1, SIZE(matrix_h, 1) 439 DO img = 1, nimg 440 CALL dbcsr_finalize(matrix_h(i, img)%matrix) 441 END DO 442 ENDDO 443 444 ! set repulsive energy 445 CALL mp_sum(erep, para_env%group) 446 energy%repulsive = erep 447 448 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 449 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 450 qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN 451 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", & 452 extension=".Log") 453 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 454 after = MIN(MAX(after, 1), 16) 455 DO img = 1, nimg 456 CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, & 457 output_unit=iw, omit_headers=omit_headers) 458 END DO 459 460 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 461 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN") 462 END IF 463 464 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 465 qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN 466 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", & 467 extension=".Log") 468 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 469 after = MIN(MAX(after, 1), 16) 470 DO img = 1, nimg 471 CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, & 472 output_unit=iw, omit_headers=omit_headers) 473 474 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 475 qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN 476 DO i = 2, SIZE(matrix_s, 1) 477 CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, & 478 output_unit=iw, omit_headers=omit_headers) 479 END DO 480 END IF 481 END DO 482 483 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 484 "DFT%PRINT%AO_MATRICES/OVERLAP") 485 END IF 486 487 IF (calculate_forces) THEN 488 IF (SIZE(matrix_p, 1) == 2) THEN 489 DO img = 1, nimg 490 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, & 491 beta_scalar=-1.0_dp) 492 END DO 493 END IF 494 DEALLOCATE (atom_of_kind) 495 END IF 496 497 CALL timestop(handle) 498 499 END SUBROUTINE build_dftb_matrices 500 501! ************************************************************************************************** 502!> \brief ... 503!> \param qs_env ... 504!> \param calculate_forces ... 505!> \param just_energy ... 506! ************************************************************************************************** 507 SUBROUTINE build_dftb_ks_matrix(qs_env, calculate_forces, just_energy) 508 TYPE(qs_environment_type), POINTER :: qs_env 509 LOGICAL, INTENT(in) :: calculate_forces, just_energy 510 511 CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb_ks_matrix', & 512 routineP = moduleN//':'//routineN 513 514 INTEGER :: atom_a, handle, iatom, ikind, img, & 515 ispin, natom, nkind, nspins, & 516 output_unit 517 REAL(KIND=dp) :: pc_ener, qmmm_el, zeff 518 REAL(KIND=dp), DIMENSION(:), POINTER :: mcharge, occupation_numbers 519 REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges 520 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 521 TYPE(cp_logger_type), POINTER :: logger 522 TYPE(cp_para_env_type), POINTER :: para_env 523 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1, mo_derivs 524 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_h, matrix_p, matrix_s 525 TYPE(dbcsr_type), POINTER :: mo_coeff 526 TYPE(dft_control_type), POINTER :: dft_control 527 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 528 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 529 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind 530 TYPE(qs_energy_type), POINTER :: energy 531 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 532 TYPE(qs_ks_env_type), POINTER :: ks_env 533 TYPE(qs_rho_type), POINTER :: rho 534 TYPE(section_vals_type), POINTER :: scf_section 535 536 CALL timeset(routineN, handle) 537 NULLIFY (dft_control, logger, scf_section, matrix_p, particle_set, ks_env, & 538 ks_matrix, rho, energy) 539 logger => cp_get_default_logger() 540 CPASSERT(ASSOCIATED(qs_env)) 541 542 CALL get_qs_env(qs_env, & 543 dft_control=dft_control, & 544 atomic_kind_set=atomic_kind_set, & 545 qs_kind_set=qs_kind_set, & 546 matrix_h_kp=matrix_h, & 547 para_env=para_env, & 548 ks_env=ks_env, & 549 matrix_ks_kp=ks_matrix, & 550 rho=rho, & 551 energy=energy) 552 553 energy%hartree = 0.0_dp 554 energy%qmmm_el = 0.0_dp 555 556 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF") 557 nspins = dft_control%nspins 558 CPASSERT(ASSOCIATED(matrix_h)) 559 CPASSERT(ASSOCIATED(rho)) 560 CPASSERT(SIZE(ks_matrix) > 0) 561 562 DO ispin = 1, nspins 563 DO img = 1, SIZE(ks_matrix, 2) 564 ! copy the core matrix into the fock matrix 565 CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix) 566 END DO 567 END DO 568 569 IF (dft_control%qs_control%dftb_control%self_consistent) THEN 570 ! Mulliken charges 571 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, & 572 matrix_s_kp=matrix_s) 573 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 574 natom = SIZE(particle_set) 575 ALLOCATE (charges(natom, nspins)) 576 ! 577 CALL mulliken_charges(matrix_p, matrix_s, para_env, charges) 578 ! 579 ALLOCATE (mcharge(natom)) 580 nkind = SIZE(atomic_kind_set) 581 DO ikind = 1, nkind 582 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom) 583 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind) 584 CALL get_dftb_atom_param(dftb_kind, zeff=zeff) 585 DO iatom = 1, natom 586 atom_a = atomic_kind_set(ikind)%atom_list(iatom) 587 mcharge(atom_a) = zeff - SUM(charges(atom_a, 1:nspins)) 588 END DO 589 END DO 590 DEALLOCATE (charges) 591 592 CALL build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, & 593 calculate_forces, just_energy) 594 595 CALL efield_tb_matrix(qs_env, ks_matrix, rho, mcharge, energy, & 596 calculate_forces, just_energy) 597 598 DEALLOCATE (mcharge) 599 600 END IF 601 602 IF (qs_env%qmmm) THEN 603 CPASSERT(SIZE(ks_matrix, 2) == 1) 604 DO ispin = 1, nspins 605 ! If QM/MM sumup the 1el Hamiltonian 606 CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, qs_env%ks_qmmm_env%matrix_h(1)%matrix, & 607 1.0_dp, 1.0_dp) 608 CALL qs_rho_get(rho, rho_ao=matrix_p1) 609 ! Compute QM/MM Energy 610 CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, & 611 matrix_p1(ispin)%matrix, qmmm_el) 612 energy%qmmm_el = energy%qmmm_el + qmmm_el 613 END DO 614 pc_ener = qs_env%ks_qmmm_env%pc_ener 615 energy%qmmm_el = energy%qmmm_el + pc_ener 616 END IF 617 618 energy%total = energy%core + energy%hartree + energy%qmmm_el + energy%efield + & 619 energy%repulsive + energy%dispersion + energy%dftb3 620 621 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", & 622 extension=".scfLog") 623 IF (output_unit > 0) THEN 624 WRITE (UNIT=output_unit, FMT="(/,(T9,A,T60,F20.10))") & 625 "Repulsive pair potential energy: ", energy%repulsive, & 626 "Zeroth order Hamiltonian energy: ", energy%core, & 627 "Charge fluctuation energy: ", energy%hartree, & 628 "London dispersion energy: ", energy%dispersion 629 IF (ABS(energy%efield) > 1.e-12_dp) THEN 630 WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") & 631 "Electric field interaction energy: ", energy%efield 632 END IF 633 IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN 634 WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") & 635 "DFTB3 3rd Order Energy Correction ", energy%dftb3 636 END IF 637 IF (qs_env%qmmm) THEN 638 WRITE (UNIT=output_unit, FMT="(T9,A,T60,F20.10)") & 639 "QM/MM Electrostatic energy: ", energy%qmmm_el 640 END IF 641 END IF 642 CALL cp_print_key_finished_output(output_unit, logger, scf_section, & 643 "PRINT%DETAILED_ENERGY") 644 ! here we compute dE/dC if needed. Assumes dE/dC is H_{ks}C (plus occupation numbers) 645 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy) THEN 646 CPASSERT(SIZE(ks_matrix, 2) == 1) 647 CALL get_qs_env(qs_env, mo_derivs=mo_derivs, mos=mo_array) 648 DO ispin = 1, SIZE(mo_derivs) 649 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, & 650 mo_coeff_b=mo_coeff, occupation_numbers=occupation_numbers) 651 IF (.NOT. mo_array(ispin)%mo_set%use_mo_coeff_b) THEN 652 CPABORT("") 653 ENDIF 654 CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin, 1)%matrix, mo_coeff, & 655 0.0_dp, mo_derivs(ispin)%matrix) 656 ENDDO 657 ENDIF 658 659 CALL timestop(handle) 660 661 END SUBROUTINE build_dftb_ks_matrix 662 663! ************************************************************************************************** 664!> \brief ... 665!> \param qs_env ... 666!> \param nderivative ... 667!> \param matrix_s ... 668! ************************************************************************************************** 669 SUBROUTINE build_dftb_overlap(qs_env, nderivative, matrix_s) 670 671 TYPE(qs_environment_type), POINTER :: qs_env 672 INTEGER, INTENT(IN) :: nderivative 673 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 674 675 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_dftb_overlap', & 676 routineP = moduleN//':'//routineN 677 678 INTEGER :: handle, i, iatom, icol, ikind, indder, irow, j, jatom, jkind, l1, l2, la, lb, & 679 llm, lmaxi, lmaxj, m, n1, n2, natom, natorb_a, natorb_b, ngrd, ngrdcut, nkind 680 LOGICAL :: defined, found 681 REAL(KIND=dp) :: ddr, dgrd, dr, f0 682 REAL(KIND=dp), DIMENSION(0:3) :: skself 683 REAL(KIND=dp), DIMENSION(3) :: drij, rij 684 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, sblock, smatij, smatji 685 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsblock1 686 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 687 TYPE(block_p_type), DIMENSION(2:10) :: dsblocks 688 TYPE(cp_logger_type), POINTER :: logger 689 TYPE(dft_control_type), POINTER :: dft_control 690 TYPE(dftb_control_type), POINTER :: dftb_control 691 TYPE(neighbor_list_iterator_p_type), & 692 DIMENSION(:), POINTER :: nl_iterator 693 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 694 POINTER :: sab_orb 695 TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b 696 TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), & 697 POINTER :: dftb_potential 698 TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji 699 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 700 701 CALL timeset(routineN, handle) 702 703 ! set pointers 704 iptr = 0 705 DO la = 0, 3 706 DO lb = 0, 3 707 llm = 0 708 DO l1 = 0, MAX(la, lb) 709 DO l2 = 0, MIN(l1, la, lb) 710 DO m = 0, l2 711 llm = llm + 1 712 iptr(l1, l2, m, la, lb) = llm 713 END DO 714 END DO 715 END DO 716 END DO 717 END DO 718 719 NULLIFY (logger) 720 logger => cp_get_default_logger() 721 722 NULLIFY (atomic_kind_set, qs_kind_set, sab_orb) 723 724 CALL get_qs_env(qs_env=qs_env, & 725 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & 726 dft_control=dft_control) 727 728 dftb_control => dft_control%qs_control%dftb_control 729 730 NULLIFY (dftb_potential) 731 CALL get_qs_env(qs_env=qs_env, & 732 dftb_potential=dftb_potential) 733 734 nkind = SIZE(atomic_kind_set) 735 736 ! Allocate the overlap matrix 737 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb) 738 CALL setup_matrices1(qs_env, nderivative, matrix_s, 'OVERLAP', sab_orb) 739 740 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 741 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 742 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 743 iatom=iatom, jatom=jatom, r=rij) 744 745 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom) 746 CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a) 747 CALL get_dftb_atom_param(dftb_kind_a, & 748 defined=defined, lmax=lmaxi, skself=skself, & 749 natorb=natorb_a) 750 751 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 752 753 CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b) 754 CALL get_dftb_atom_param(dftb_kind_b, & 755 defined=defined, lmax=lmaxj, natorb=natorb_b) 756 757 IF (.NOT. defined .OR. natorb_b < 1) CYCLE 758 759 ! retrieve information on F and S matrix 760 dftb_param_ij => dftb_potential(ikind, jkind) 761 dftb_param_ji => dftb_potential(jkind, ikind) 762 ! assume table size and type is symmetric 763 ngrd = dftb_param_ij%ngrd 764 ngrdcut = dftb_param_ij%ngrdcut 765 dgrd = dftb_param_ij%dgrd 766 ddr = dgrd*0.1_dp 767 CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm) 768 llm = dftb_param_ij%llm 769 smatij => dftb_param_ij%smat 770 smatji => dftb_param_ji%smat 771 772 dr = SQRT(SUM(rij(:)**2)) 773 IF (NINT(dr/dgrd) <= ngrdcut) THEN 774 775 icol = MAX(iatom, jatom); irow = MIN(iatom, jatom) 776 777 NULLIFY (sblock) 778 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, & 779 row=irow, col=icol, BLOCK=sblock, found=found) 780 CPASSERT(found) 781 782 IF (nderivative .GT. 0) THEN 783 DO i = 2, SIZE(matrix_s, 1) 784 NULLIFY (dsblocks(i)%block) 785 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, & 786 row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found) 787 END DO 788 END IF 789 790 IF (iatom == jatom .AND. dr < 0.001_dp) THEN 791 ! diagonal block 792 DO i = 1, natorb_a 793 sblock(i, i) = sblock(i, i) + 1._dp 794 END DO 795 ELSE 796 ! off-diagonal block 797 CALL compute_block_sk(sblock, smatij, smatji, rij, ngrd, ngrdcut, dgrd, & 798 llm, lmaxi, lmaxj, irow, iatom) 799 800 IF (nderivative .GE. 1) THEN 801 n1 = SIZE(sblock, 1); n2 = SIZE(sblock, 2) 802 indder = 1 ! used to put the 2nd derivatives in the correct matric (5=xx,8=yy,10=zz) 803 804 ALLOCATE (dsblock1(n1, n2, 3), dsblock(n1, n2), dsblockm(n1, n2)) 805 dsblock1 = 0.0_dp 806 DO i = 1, 3 807 dsblock = 0._dp; dsblockm = 0.0_dp 808 drij = rij 809 f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp 810 811 drij(i) = rij(i) - ddr*f0 812 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 813 llm, lmaxi, lmaxj, irow, iatom) 814 815 drij(i) = rij(i) + ddr*f0 816 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 817 llm, lmaxi, lmaxj, irow, iatom) 818 819 dsblock1(:, :, i) = (dsblock + dsblockm) 820 dsblock = dsblock - dsblockm 821 dsblock = dsblock/(2.0_dp*ddr) 822 823 CPASSERT(ASSOCIATED(dsblocks(i + 1)%block)) 824 dsblocks(i + 1)%block = dsblocks(i + 1)%block + dsblock 825 IF (nderivative .GT. 1) THEN 826 indder = indder + 5 - i 827 dsblocks(indder)%block = 0.0_dp 828 dsblocks(indder)%block = dsblocks(indder)%block + & 829 (dsblock1(:, :, i) - 2.0_dp*sblock)/ddr**2 830 END IF 831 ENDDO 832 833 IF (nderivative .GT. 1) THEN 834 DO i = 1, 2 835 DO j = i + 1, 3 836 dsblock = 0._dp; dsblockm = 0.0_dp 837 drij = rij 838 f0 = 1.0_dp; IF (irow == iatom) f0 = -1.0_dp 839 840 drij(i) = rij(i) - ddr*f0; drij(j) = rij(j) - ddr*f0 841 CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 842 llm, lmaxi, lmaxj, irow, iatom) 843 844 drij(i) = rij(i) + ddr*f0; drij(j) = rij(j) + ddr*f0 845 CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, & 846 llm, lmaxi, lmaxj, irow, iatom) 847 848 indder = 2 + 2*i + j 849 dsblocks(indder)%block = 0.0_dp 850 dsblocks(indder)%block = & 851 dsblocks(indder)%block + ( & 852 dsblock + dsblockm - dsblock1(:, :, i) - dsblock1(:, :, j) + 2.0_dp*sblock)/(2.0_dp*ddr**2) 853 END DO 854 END DO 855 END IF 856 857 DEALLOCATE (dsblock1, dsblock, dsblockm) 858 END IF 859 END IF 860 END IF 861 END DO 862 CALL neighbor_list_iterator_release(nl_iterator) 863 864 DO i = 1, SIZE(matrix_s, 1) 865 CALL dbcsr_finalize(matrix_s(i)%matrix) 866 ENDDO 867 868 CALL timestop(handle) 869 870 END SUBROUTINE build_dftb_overlap 871 872! ************************************************************************************************** 873!> \brief ... 874!> \param qs_env ... 875!> \param nderivative ... 876!> \param matrices ... 877!> \param mnames ... 878!> \param sab_nl ... 879! ************************************************************************************************** 880 SUBROUTINE setup_matrices1(qs_env, nderivative, matrices, mnames, sab_nl) 881 882 TYPE(qs_environment_type), POINTER :: qs_env 883 INTEGER, INTENT(IN) :: nderivative 884 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices 885 CHARACTER(LEN=*) :: mnames 886 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 887 POINTER :: sab_nl 888 889 CHARACTER(len=*), PARAMETER :: routineN = 'setup_matrices1', & 890 routineP = moduleN//':'//routineN 891 892 CHARACTER(1) :: symmetry_type 893 CHARACTER(LEN=default_string_length) :: matnames 894 INTEGER :: i, natom, neighbor_list_id, nkind, nmat, & 895 nsgf 896 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 897 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes 898 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 899 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 900 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 901 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 902 903 NULLIFY (particle_set, atomic_kind_set) 904 905 CALL get_qs_env(qs_env=qs_env, & 906 atomic_kind_set=atomic_kind_set, & 907 qs_kind_set=qs_kind_set, & 908 particle_set=particle_set, & 909 dbcsr_dist=dbcsr_dist, & 910 neighbor_list_id=neighbor_list_id) 911 912 nkind = SIZE(atomic_kind_set) 913 natom = SIZE(particle_set) 914 915 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf) 916 917 ALLOCATE (first_sgf(natom)) 918 ALLOCATE (last_sgf(natom)) 919 920 CALL get_particle_set(particle_set, qs_kind_set, & 921 first_sgf=first_sgf, & 922 last_sgf=last_sgf) 923 924 nmat = 0 925 IF (nderivative == 0) nmat = 1 926 IF (nderivative == 1) nmat = 4 927 IF (nderivative == 2) nmat = 10 928 CPASSERT(nmat > 0) 929 930 ALLOCATE (row_blk_sizes(natom)) 931 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf) 932 933 CALL dbcsr_allocate_matrix_set(matrices, nmat) 934 935 ! Up to 2nd derivative take care to get the symmetries correct 936 DO i = 1, nmat 937 IF (i .GT. 1) THEN 938 matnames = TRIM(mnames)//" DERIVATIVE MATRIX DFTB" 939 symmetry_type = dbcsr_type_antisymmetric 940 IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric 941 ELSE 942 symmetry_type = dbcsr_type_symmetric 943 matnames = TRIM(mnames)//" MATRIX DFTB" 944 END IF 945 ALLOCATE (matrices(i)%matrix) 946 CALL dbcsr_create(matrix=matrices(i)%matrix, & 947 name=TRIM(matnames), & 948 dist=dbcsr_dist, matrix_type=symmetry_type, & 949 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 950 nze=0, mutable_work=.TRUE.) 951 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i)%matrix, sab_nl) 952 END DO 953 954 DEALLOCATE (first_sgf) 955 DEALLOCATE (last_sgf) 956 957 DEALLOCATE (row_blk_sizes) 958 959 END SUBROUTINE setup_matrices1 960 961! ************************************************************************************************** 962!> \brief ... 963!> \param qs_env ... 964!> \param nderivative ... 965!> \param nimg ... 966!> \param matrices ... 967!> \param mnames ... 968!> \param sab_nl ... 969! ************************************************************************************************** 970 SUBROUTINE setup_matrices2(qs_env, nderivative, nimg, matrices, mnames, sab_nl) 971 972 TYPE(qs_environment_type), POINTER :: qs_env 973 INTEGER, INTENT(IN) :: nderivative, nimg 974 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrices 975 CHARACTER(LEN=*) :: mnames 976 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 977 POINTER :: sab_nl 978 979 CHARACTER(len=*), PARAMETER :: routineN = 'setup_matrices2', & 980 routineP = moduleN//':'//routineN 981 982 CHARACTER(1) :: symmetry_type 983 CHARACTER(LEN=default_string_length) :: matnames 984 INTEGER :: i, img, natom, neighbor_list_id, nkind, & 985 nmat, nsgf 986 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf 987 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes 988 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 989 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 990 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 991 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 992 993 NULLIFY (particle_set, atomic_kind_set) 994 995 CALL get_qs_env(qs_env=qs_env, & 996 atomic_kind_set=atomic_kind_set, & 997 qs_kind_set=qs_kind_set, & 998 particle_set=particle_set, & 999 dbcsr_dist=dbcsr_dist, & 1000 neighbor_list_id=neighbor_list_id) 1001 1002 nkind = SIZE(atomic_kind_set) 1003 natom = SIZE(particle_set) 1004 1005 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf) 1006 1007 ALLOCATE (first_sgf(natom)) 1008 ALLOCATE (last_sgf(natom)) 1009 1010 CALL get_particle_set(particle_set, qs_kind_set, & 1011 first_sgf=first_sgf, & 1012 last_sgf=last_sgf) 1013 1014 nmat = 0 1015 IF (nderivative == 0) nmat = 1 1016 IF (nderivative == 1) nmat = 4 1017 IF (nderivative == 2) nmat = 10 1018 CPASSERT(nmat > 0) 1019 1020 ALLOCATE (row_blk_sizes(natom)) 1021 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf) 1022 1023 CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg) 1024 1025 ! Up to 2nd derivative take care to get the symmetries correct 1026 DO img = 1, nimg 1027 DO i = 1, nmat 1028 IF (i .GT. 1) THEN 1029 matnames = TRIM(mnames)//" DERIVATIVE MATRIX DFTB" 1030 symmetry_type = dbcsr_type_antisymmetric 1031 IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric 1032 ELSE 1033 symmetry_type = dbcsr_type_symmetric 1034 matnames = TRIM(mnames)//" MATRIX DFTB" 1035 END IF 1036 ALLOCATE (matrices(i, img)%matrix) 1037 CALL dbcsr_create(matrix=matrices(i, img)%matrix, & 1038 name=TRIM(matnames), & 1039 dist=dbcsr_dist, matrix_type=symmetry_type, & 1040 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 1041 nze=0, mutable_work=.TRUE.) 1042 CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl) 1043 END DO 1044 END DO 1045 1046 DEALLOCATE (first_sgf) 1047 DEALLOCATE (last_sgf) 1048 1049 DEALLOCATE (row_blk_sizes) 1050 1051 END SUBROUTINE setup_matrices2 1052 1053END MODULE qs_dftb_matrices 1054 1055