1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Interface between ALMO SCF and QS 8!> \par History 9!> 2011.05 created [Rustam Z Khaliullin] 10!> \author Rustam Z Khaliullin 11! ************************************************************************************************** 12MODULE almo_scf_qs 13 USE almo_scf_types, ONLY: almo_mat_dim_aobasis,& 14 almo_mat_dim_occ,& 15 almo_mat_dim_virt,& 16 almo_mat_dim_virt_disc,& 17 almo_mat_dim_virt_full,& 18 almo_scf_env_type 19 USE atomic_kind_types, ONLY: get_atomic_kind 20 USE cell_types, ONLY: cell_type,& 21 pbc 22 USE cp_control_types, ONLY: dft_control_type 23 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 24 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 25 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 26 cp_fm_struct_release,& 27 cp_fm_struct_type 28 USE cp_fm_types, ONLY: cp_fm_create,& 29 cp_fm_release,& 30 cp_fm_type 31 USE cp_log_handling, ONLY: cp_get_default_logger,& 32 cp_logger_get_default_unit_nr,& 33 cp_logger_type 34 USE cp_units, ONLY: cp_unit_to_cp2k 35 USE dbcsr_api, ONLY: & 36 dbcsr_complete_redistribute, dbcsr_copy, dbcsr_copy_into_existing, dbcsr_create, & 37 dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_new, & 38 dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, & 39 dbcsr_get_block_p, dbcsr_get_info, dbcsr_get_num_blocks, dbcsr_get_stored_coordinates, & 40 dbcsr_multiply, dbcsr_nblkcols_total, dbcsr_nblkrows_total, dbcsr_p_type, dbcsr_release, & 41 dbcsr_reserve_block2d, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_work_create 42 USE input_constants, ONLY: almo_constraint_ao_overlap,& 43 almo_constraint_block_diagonal,& 44 almo_constraint_distance,& 45 almo_domain_layout_molecular,& 46 almo_mat_distr_atomic,& 47 almo_mat_distr_molecular,& 48 do_bondparm_covalent,& 49 do_bondparm_vdw 50 USE kinds, ONLY: dp 51 USE message_passing, ONLY: mp_allgather 52 USE molecule_types, ONLY: get_molecule_set_info,& 53 molecule_type 54 USE particle_types, ONLY: particle_type 55 USE qs_energy_types, ONLY: qs_energy_type 56 USE qs_environment_types, ONLY: get_qs_env,& 57 qs_environment_type,& 58 set_qs_env 59 USE qs_ks_methods, ONLY: qs_ks_update_qs_env 60 USE qs_ks_types, ONLY: qs_ks_did_change,& 61 qs_ks_env_type,& 62 set_ks_env 63 USE qs_mo_types, ONLY: allocate_mo_set,& 64 init_mo_set,& 65 mo_set_p_type 66 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 67 neighbor_list_iterate,& 68 neighbor_list_iterator_create,& 69 neighbor_list_iterator_p_type,& 70 neighbor_list_iterator_release,& 71 neighbor_list_set_p_type 72 USE qs_rho_methods, ONLY: qs_rho_update_rho 73 USE qs_rho_types, ONLY: qs_rho_get,& 74 qs_rho_type 75 USE qs_scf_types, ONLY: qs_scf_env_type,& 76 scf_env_create 77#include "./base/base_uses.f90" 78 79 IMPLICIT NONE 80 81 PRIVATE 82 83 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_qs' 84 85 PUBLIC :: matrix_almo_create, & 86 almo_scf_construct_quencher, & 87 calculate_w_matrix_almo, & 88 init_almo_ks_matrix_via_qs, & 89 almo_scf_update_ks_energy, & 90 construct_qs_mos, & 91 matrix_qs_to_almo, & 92 almo_dm_to_almo_ks, & 93 almo_dm_to_qs_env 94 95CONTAINS 96 97! ************************************************************************************************** 98!> \brief create the ALMO matrix templates 99!> \param matrix_new ... 100!> \param matrix_qs ... 101!> \param almo_scf_env ... 102!> \param name_new ... 103!> \param size_keys ... 104!> \param symmetry_new ... 105!> \param spin_key ... 106!> \param init_domains ... 107!> \par History 108!> 2011.05 created [Rustam Z Khaliullin] 109!> \author Rustam Z Khaliullin 110! ************************************************************************************************** 111 SUBROUTINE matrix_almo_create(matrix_new, matrix_qs, almo_scf_env, & 112 name_new, size_keys, symmetry_new, & 113 spin_key, init_domains) 114 115 TYPE(dbcsr_type) :: matrix_new, matrix_qs 116 TYPE(almo_scf_env_type), INTENT(IN) :: almo_scf_env 117 CHARACTER(len=*), INTENT(IN) :: name_new 118 INTEGER, DIMENSION(2), INTENT(IN) :: size_keys 119 CHARACTER, INTENT(IN) :: symmetry_new 120 INTEGER, INTENT(IN) :: spin_key 121 LOGICAL, INTENT(IN) :: init_domains 122 123 CHARACTER(len=*), PARAMETER :: routineN = 'matrix_almo_create', & 124 routineP = moduleN//':'//routineN 125 126 INTEGER :: dimen, handle, hold, iatom, iblock_col, & 127 iblock_row, imol, mynode, natoms, & 128 nblkrows_tot, nlength, nmols, row 129 INTEGER, DIMENSION(:), POINTER :: blk_distr, blk_sizes, block_sizes_new, col_distr_new, & 130 col_sizes_new, distr_new_array, row_distr_new, row_sizes_new 131 LOGICAL :: active, one_dim_is_mo, tr 132 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_new_block 133 TYPE(dbcsr_distribution_type) :: dist_new, dist_qs 134 135! dimension size: AO, MO, etc 136! almo_mat_dim_aobasis - no. of AOs, 137! almo_mat_dim_occ - no. of occupied MOs 138! almo_mat_dim_domains - no. of domains 139! symmetry type: dbcsr_type_no_symmetry, dbcsr_type_symmetric, 140! dbcsr_type_antisymmetric, dbcsr_type_hermitian, dbcsr_type_antihermitian 141! (see dbcsr_lib/dbcsr_types.F for other values) 142! spin_key: either 1 or 2 (0 is allowed for matrics in the AO basis) 143! TYPE(dbcsr_iterator_type) :: iter 144! REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: allones 145!----------------------------------------------------------------------- 146 147 CALL timeset(routineN, handle) 148 149 ! RZK-warning The structure of the matrices can be optimized: 150 ! 1. Diagonal matrices must be distributed evenly over the processes. 151 ! This can be achieved by distributing cpus: 012012-rows and 001122-cols 152 ! block_diagonal_flag is introduced but not used 153 ! 2. Multiplication of diagonally dominant matrices will be faster 154 ! if the diagonal blocks are local to the same processes. 155 ! 3. Systems of molecules of drastically different sizes might need 156 ! better distribution. 157 158 ! obtain distribution from the qs matrix - it might be useful 159 ! to get the structure of the AO dimensions 160 CALL dbcsr_get_info(matrix_qs, distribution=dist_qs) 161 162 natoms = almo_scf_env%natoms 163 nmols = almo_scf_env%nmolecules 164 165 DO dimen = 1, 2 ! 1 - row, 2 - column dimension 166 167 ! distribution pattern is the same for all matrix types (ao, occ, virt) 168 IF (dimen == 1) THEN !rows 169 CALL dbcsr_distribution_get(dist_qs, row_dist=blk_distr) 170 ELSE !columns 171 CALL dbcsr_distribution_get(dist_qs, col_dist=blk_distr) 172 ENDIF 173 174 IF (size_keys(dimen) == almo_mat_dim_aobasis) THEN ! this dimension is AO 175 176 ! structure of an AO dimension can be copied from matrix_qs 177 CALL dbcsr_get_info(matrix_qs, row_blk_size=blk_sizes) 178 179 ! atomic clustering of AOs 180 IF (almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN 181 ALLOCATE (block_sizes_new(natoms), distr_new_array(natoms)) 182 block_sizes_new(:) = blk_sizes(:) 183 distr_new_array(:) = blk_distr(:) 184 ! molecular clustering of AOs 185 ELSE IF (almo_scf_env%mat_distr_aos == almo_mat_distr_molecular) THEN 186 ALLOCATE (block_sizes_new(nmols), distr_new_array(nmols)) 187 block_sizes_new(:) = 0 188 DO iatom = 1, natoms 189 block_sizes_new(almo_scf_env%domain_index_of_atom(iatom)) = & 190 block_sizes_new(almo_scf_env%domain_index_of_atom(iatom)) + & 191 blk_sizes(iatom) 192 ENDDO 193 DO imol = 1, nmols 194 distr_new_array(imol) = & 195 blk_distr(almo_scf_env%first_atom_of_domain(imol)) 196 ENDDO 197 ELSE 198 CPABORT("Illegal distribution") 199 ENDIF 200 201 ELSE ! this dimension is not AO 202 203 IF (size_keys(dimen) == almo_mat_dim_occ .OR. & 204 size_keys(dimen) == almo_mat_dim_virt .OR. & 205 size_keys(dimen) == almo_mat_dim_virt_disc .OR. & 206 size_keys(dimen) == almo_mat_dim_virt_full) THEN ! this dim is MO 207 208 ! atomic clustering of MOs 209 IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN 210 nlength = natoms 211 ALLOCATE (block_sizes_new(nlength)) 212 block_sizes_new(:) = 0 213 IF (size_keys(dimen) == almo_mat_dim_occ) THEN 214 ! currently distributing atomic distr of mos is not allowed 215 ! RZK-warning define nocc_of_atom and nvirt_atom to implement it 216 !block_sizes_new(:)=almo_scf_env%nocc_of_atom(:,spin_key) 217 ELSE IF (size_keys(dimen) == almo_mat_dim_virt) THEN 218 !block_sizes_new(:)=almo_scf_env%nvirt_of_atom(:,spin_key) 219 ENDIF 220 ! molecular clustering of MOs 221 ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN 222 nlength = nmols 223 ALLOCATE (block_sizes_new(nlength)) 224 IF (size_keys(dimen) == almo_mat_dim_occ) THEN 225 block_sizes_new(:) = almo_scf_env%nocc_of_domain(:, spin_key) 226 ! Handle zero-electron fragments by adding one-orbital that 227 ! must remain zero at all times 228 WHERE (block_sizes_new == 0) block_sizes_new = 1 229 ELSE IF (size_keys(dimen) == almo_mat_dim_virt_disc) THEN 230 block_sizes_new(:) = almo_scf_env%nvirt_disc_of_domain(:, spin_key) 231 ELSE IF (size_keys(dimen) == almo_mat_dim_virt_full) THEN 232 block_sizes_new(:) = almo_scf_env%nvirt_full_of_domain(:, spin_key) 233 ELSE IF (size_keys(dimen) == almo_mat_dim_virt) THEN 234 block_sizes_new(:) = almo_scf_env%nvirt_of_domain(:, spin_key) 235 ENDIF 236 ELSE 237 CPABORT("Illegal distribution") 238 ENDIF 239 240 ELSE 241 242 CPABORT("Illegal dimension") 243 244 ENDIF ! end choosing dim size (occ, virt) 245 246 ! distribution for MOs is copied from AOs 247 ALLOCATE (distr_new_array(nlength)) 248 ! atomic clustering 249 IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN 250 distr_new_array(:) = blk_distr(:) 251 ! molecular clustering 252 ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN 253 DO imol = 1, nmols 254 distr_new_array(imol) = & 255 blk_distr(almo_scf_env%first_atom_of_domain(imol)) 256 ENDDO 257 ENDIF 258 ENDIF ! end choosing dimension size (AOs vs .NOT.AOs) 259 260 ! create final arrays 261 IF (dimen == 1) THEN !rows 262 row_sizes_new => block_sizes_new 263 row_distr_new => distr_new_array 264 ELSE !columns 265 col_sizes_new => block_sizes_new 266 col_distr_new => distr_new_array 267 ENDIF 268 ENDDO ! both rows and columns are done 269 270 ! Create the distribution 271 CALL dbcsr_distribution_new(dist_new, template=dist_qs, & 272 row_dist=row_distr_new, col_dist=col_distr_new, & 273 reuse_arrays=.TRUE.) 274 275 ! Create the matrix 276 CALL dbcsr_create(matrix_new, name_new, & 277 dist_new, symmetry_new, & 278 row_sizes_new, col_sizes_new, reuse_arrays=.TRUE.) 279 CALL dbcsr_distribution_release(dist_new) 280 281 ! fill out reqired blocks with 1.0_dp to tell the dbcsr library 282 ! which blocks to keep 283 IF (init_domains) THEN 284 285 CALL dbcsr_distribution_get(dist_new, mynode=mynode) 286 CALL dbcsr_work_create(matrix_new, work_mutable=.TRUE.) 287 288 ! startQQQ - this part of the code scales quadratically 289 ! therefore it is replaced with a less general but linear scaling algorithm below 290 ! the quadratic algorithm is kept to be re-written later 291 !QQQnblkrows_tot = dbcsr_nblkrows_total(matrix_new) 292 !QQQnblkcols_tot = dbcsr_nblkcols_total(matrix_new) 293 !QQQDO row = 1, nblkrows_tot 294 !QQQ DO col = 1, nblkcols_tot 295 !QQQ tr = .FALSE. 296 !QQQ iblock_row = row 297 !QQQ iblock_col = col 298 !QQQ CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, tr, hold) 299 300 !QQQ IF(hold.EQ.mynode) THEN 301 !QQQ 302 !QQQ ! RZK-warning replace with a function which says if this 303 !QQQ ! distribution block is active or not 304 !QQQ ! Translate indeces of distribution blocks to domain blocks 305 !QQQ if (size_keys(1)==almo_mat_dim_aobasis) then 306 !QQQ domain_row=almo_scf_env%domain_index_of_ao_block(iblock_row) 307 !QQQ else if (size_keys(2)==almo_mat_dim_occ .OR. & 308 !QQQ size_keys(2)==almo_mat_dim_virt .OR. & 309 !QQQ size_keys(2)==almo_mat_dim_virt_disc .OR. & 310 !QQQ size_keys(2)==almo_mat_dim_virt_full) then 311 !QQQ domain_row=almo_scf_env%domain_index_of_mo_block(iblock_row) 312 !QQQ else 313 !QQQ CPErrorMessage(cp_failure_level,routineP,"Illegal dimension") 314 !QQQ CPPrecondition(.FALSE.,cp_failure_level,routineP,failure) 315 !QQQ endif 316 317 !QQQ if (size_keys(2)==almo_mat_dim_aobasis) then 318 !QQQ domain_col=almo_scf_env%domain_index_of_ao_block(iblock_col) 319 !QQQ else if (size_keys(2)==almo_mat_dim_occ .OR. & 320 !QQQ size_keys(2)==almo_mat_dim_virt .OR. & 321 !QQQ size_keys(2)==almo_mat_dim_virt_disc .OR. & 322 !QQQ size_keys(2)==almo_mat_dim_virt_full) then 323 !QQQ domain_col=almo_scf_env%domain_index_of_mo_block(iblock_col) 324 !QQQ else 325 !QQQ CPErrorMessage(cp_failure_level,routineP,"Illegal dimension") 326 !QQQ CPPrecondition(.FALSE.,cp_failure_level,routineP,failure) 327 !QQQ endif 328 329 !QQQ ! Finds if we need this block 330 !QQQ ! only the block-diagonal constraint is implemented here 331 !QQQ active=.false. 332 !QQQ if (domain_row==domain_col) active=.true. 333 334 !QQQ IF (active) THEN 335 !QQQ NULLIFY (p_new_block) 336 !QQQ CALL dbcsr_reserve_block2d(matrix_new, iblock_row, iblock_col, p_new_block) 337 !QQQ CPPostcondition(ASSOCIATED(p_new_block),cp_failure_level,routineP,failure) 338 !QQQ p_new_block(:,:) = 1.0_dp 339 !QQQ ENDIF 340 341 !QQQ ENDIF ! mynode 342 !QQQ ENDDO 343 !QQQENDDO 344 !QQQtake care of zero-electron fragments 345 ! endQQQ - end of the quadratic part 346 ! start linear-scaling replacement: 347 ! works only for molecular blocks AND molecular distributions 348 nblkrows_tot = dbcsr_nblkrows_total(matrix_new) 349 DO row = 1, nblkrows_tot 350 tr = .FALSE. 351 iblock_row = row 352 iblock_col = row 353 CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, hold) 354 355 IF (hold .EQ. mynode) THEN 356 357 active = .TRUE. 358 359 one_dim_is_mo = .FALSE. 360 DO dimen = 1, 2 ! 1 - row, 2 - column dimension 361 IF (size_keys(dimen) == almo_mat_dim_occ) one_dim_is_mo = .TRUE. 362 ENDDO 363 IF (one_dim_is_mo) THEN 364 IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .FALSE. 365 END IF 366 367 one_dim_is_mo = .FALSE. 368 DO dimen = 1, 2 369 IF (size_keys(dimen) == almo_mat_dim_virt) one_dim_is_mo = .TRUE. 370 ENDDO 371 IF (one_dim_is_mo) THEN 372 IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .FALSE. 373 END IF 374 375 one_dim_is_mo = .FALSE. 376 DO dimen = 1, 2 377 IF (size_keys(dimen) == almo_mat_dim_virt_disc) one_dim_is_mo = .TRUE. 378 ENDDO 379 IF (one_dim_is_mo) THEN 380 IF (almo_scf_env%nvirt_disc_of_domain(row, spin_key) == 0) active = .FALSE. 381 END IF 382 383 one_dim_is_mo = .FALSE. 384 DO dimen = 1, 2 385 IF (size_keys(dimen) == almo_mat_dim_virt_full) one_dim_is_mo = .TRUE. 386 ENDDO 387 IF (one_dim_is_mo) THEN 388 IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .FALSE. 389 END IF 390 391 IF (active) THEN 392 NULLIFY (p_new_block) 393 CALL dbcsr_reserve_block2d(matrix_new, iblock_row, iblock_col, p_new_block) 394 CPASSERT(ASSOCIATED(p_new_block)) 395 p_new_block(:, :) = 1.0_dp 396 ENDIF 397 398 ENDIF ! mynode 399 ENDDO 400 ! end lnear-scaling replacement 401 402 ENDIF ! init_domains 403 404 CALL dbcsr_finalize(matrix_new) 405 406 CALL timestop(handle) 407 408 END SUBROUTINE matrix_almo_create 409 410! ************************************************************************************************** 411!> \brief convert between two types of matrices: QS style to ALMO style 412!> \param matrix_qs ... 413!> \param matrix_almo ... 414!> \param mat_distr_aos ... 415!> \param keep_sparsity ... 416!> \par History 417!> 2011.06 created [Rustam Z Khaliullin] 418!> \author Rustam Z Khaliullin 419! ************************************************************************************************** 420 SUBROUTINE matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos, keep_sparsity) 421 422 TYPE(dbcsr_type) :: matrix_qs, matrix_almo 423 INTEGER :: mat_distr_aos 424 LOGICAL, INTENT(IN) :: keep_sparsity 425 426 CHARACTER(len=*), PARAMETER :: routineN = 'matrix_qs_to_almo', & 427 routineP = moduleN//':'//routineN 428 429 INTEGER :: handle 430 TYPE(dbcsr_type) :: matrix_qs_nosym 431 432 CALL timeset(routineN, handle) 433 !RZK-warning if it's not a N(AO)xN(AO) matrix then stop 434 435 SELECT CASE (mat_distr_aos) 436 CASE (almo_mat_distr_atomic) 437 ! automatic data_type conversion 438 CALL dbcsr_copy(matrix_almo, matrix_qs, & 439 keep_sparsity=keep_sparsity) 440 CASE (almo_mat_distr_molecular) 441 ! desymmetrize the qs matrix 442 CALL dbcsr_create(matrix_qs_nosym, template=matrix_qs, & 443 matrix_type=dbcsr_type_no_symmetry) 444 CALL dbcsr_desymmetrize(matrix_qs, matrix_qs_nosym) 445 446 ! perform the magic complete_redistribute 447 ! before calling complete_redistribute set all blocks to zero 448 ! otherwise the non-zero elements of the redistributed matrix, 449 ! which are in zero-blocks of the original matrix, will remain 450 ! in the final redistributed matrix. this is a bug in 451 ! complete_redistribute. RZK-warning it should be later corrected by calling 452 ! dbcsr_set to 0.0 from within complete_redistribute 453 CALL dbcsr_set(matrix_almo, 0.0_dp) 454 CALL dbcsr_complete_redistribute(matrix_qs_nosym, matrix_almo, & 455 keep_sparsity=keep_sparsity); 456 CALL dbcsr_release(matrix_qs_nosym) 457 458 CASE DEFAULT 459 CPABORT("") 460 END SELECT 461 462 CALL timestop(handle) 463 464 END SUBROUTINE matrix_qs_to_almo 465 466! ************************************************************************************************** 467!> \brief convert between two types of matrices: ALMO style to QS style 468!> \param matrix_almo ... 469!> \param matrix_qs ... 470!> \param mat_distr_aos ... 471!> \par History 472!> 2011.06 created [Rustam Z Khaliullin] 473!> \author Rustam Z Khaliullin 474! ************************************************************************************************** 475 SUBROUTINE matrix_almo_to_qs(matrix_almo, matrix_qs, mat_distr_aos) 476 TYPE(dbcsr_type) :: matrix_almo, matrix_qs 477 INTEGER, INTENT(IN) :: mat_distr_aos 478 479 CHARACTER(len=*), PARAMETER :: routineN = 'matrix_almo_to_qs', & 480 routineP = moduleN//':'//routineN 481 482 INTEGER :: handle 483 484 CALL timeset(routineN, handle) 485 ! RZK-warning if it's not a N(AO)xN(AO) matrix then stop 486 487 SELECT CASE (mat_distr_aos) 488 CASE (almo_mat_distr_atomic) 489 CALL dbcsr_copy_into_existing(matrix_qs, matrix_almo) 490 CASE (almo_mat_distr_molecular) 491 CALL dbcsr_set(matrix_qs, 0.0_dp) 492 CALL dbcsr_complete_redistribute(matrix_almo, matrix_qs, keep_sparsity=.TRUE.) 493 CASE DEFAULT 494 CPABORT("") 495 END SELECT 496 497 CALL timestop(handle) 498 499 END SUBROUTINE matrix_almo_to_qs 500 501! ************************************************************************************************** 502!> \brief Initialization of the QS and ALMO KS matrix 503!> \param qs_env ... 504!> \param matrix_ks ... 505!> \param mat_distr_aos ... 506!> \param eps_filter ... 507!> \par History 508!> 2011.05 created [Rustam Z Khaliullin] 509!> \author Rustam Z Khaliullin 510! ************************************************************************************************** 511 SUBROUTINE init_almo_ks_matrix_via_qs(qs_env, matrix_ks, mat_distr_aos, eps_filter) 512 513 TYPE(qs_environment_type), POINTER :: qs_env 514 TYPE(dbcsr_type), DIMENSION(:) :: matrix_ks 515 INTEGER :: mat_distr_aos 516 REAL(KIND=dp) :: eps_filter 517 518 CHARACTER(len=*), PARAMETER :: routineN = 'init_almo_ks_matrix_via_qs', & 519 routineP = moduleN//':'//routineN 520 521 INTEGER :: handle, ispin, nspin 522 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks, matrix_qs_s 523 TYPE(dft_control_type), POINTER :: dft_control 524 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 525 POINTER :: sab_orb 526 TYPE(qs_ks_env_type), POINTER :: ks_env 527 528 CALL timeset(routineN, handle) 529 530 NULLIFY (sab_orb) 531 532 ! get basic quantities from the qs_env 533 CALL get_qs_env(qs_env, & 534 dft_control=dft_control, & 535 matrix_s=matrix_qs_s, & 536 matrix_ks=matrix_qs_ks, & 537 ks_env=ks_env, & 538 sab_orb=sab_orb) 539 540 nspin = dft_control%nspins 541 542 ! create matrix_ks in the QS env if necessary 543 IF (.NOT. ASSOCIATED(matrix_qs_ks)) THEN 544 CALL dbcsr_allocate_matrix_set(matrix_qs_ks, nspin) 545 DO ispin = 1, nspin 546 ALLOCATE (matrix_qs_ks(ispin)%matrix) 547 CALL dbcsr_create(matrix_qs_ks(ispin)%matrix, & 548 template=matrix_qs_s(1)%matrix) 549 CALL cp_dbcsr_alloc_block_from_nbl(matrix_qs_ks(ispin)%matrix, sab_orb) 550 CALL dbcsr_set(matrix_qs_ks(ispin)%matrix, 0.0_dp) 551 ENDDO 552 CALL set_ks_env(ks_env, matrix_ks=matrix_qs_ks) 553 ENDIF 554 555 ! copy to ALMO 556 DO ispin = 1, nspin 557 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, & 558 matrix_ks(ispin), mat_distr_aos, .FALSE.) 559 CALL dbcsr_filter(matrix_ks(ispin), eps_filter) 560 ENDDO 561 562 CALL timestop(handle) 563 564 END SUBROUTINE init_almo_ks_matrix_via_qs 565 566! ************************************************************************************************** 567!> \brief Create MOs in the QS env to be able to return ALMOs to QS 568!> \param qs_env ... 569!> \param almo_scf_env ... 570!> \par History 571!> 2016.12 created [Yifei Shi] 572!> \author Yifei Shi 573! ************************************************************************************************** 574 SUBROUTINE construct_qs_mos(qs_env, almo_scf_env) 575 576 TYPE(qs_environment_type), POINTER :: qs_env 577 TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env 578 579 CHARACTER(len=*), PARAMETER :: routineN = 'construct_qs_mos', & 580 routineP = moduleN//':'//routineN 581 582 INTEGER :: handle, ispin, ncol_fm, nrow_fm 583 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 584 TYPE(cp_fm_type), POINTER :: mo_fm_copy 585 TYPE(dft_control_type), POINTER :: dft_control 586 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 587 TYPE(qs_scf_env_type), POINTER :: scf_env 588 589 CALL timeset(routineN, handle) 590 591 ! create and init scf_env (this is necessary to return MOs to qs) 592 NULLIFY (mos, mo_fm_copy, fm_struct_tmp, scf_env) 593 CALL scf_env_create(scf_env) 594 595 !CALL qs_scf_env_initialize(qs_env, scf_env) 596 CALL set_qs_env(qs_env, scf_env=scf_env) 597 CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos) 598 599 CALL dbcsr_get_info(almo_scf_env%matrix_t(1), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm) 600 601 ! allocate and init mo_set 602 DO ispin = 1, almo_scf_env%nspins 603 604 ! Currently only fm version of mo_set is usable. 605 ! First transform the matrix_t to fm version 606 CALL allocate_mo_set(mo_set=mos(ispin)%mo_set, & 607 nao=nrow_fm, & 608 nmo=ncol_fm, & 609 nelectron=almo_scf_env%nelectrons_total, & 610 n_el_f=REAL(almo_scf_env%nelectrons_total, dp), & 611 maxocc=2.0_dp, & 612 flexible_electron_count=dft_control%relax_multiplicity) 613 614 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_fm, ncol_global=ncol_fm, & 615 context=almo_scf_env%blacs_env, & 616 para_env=almo_scf_env%para_env) 617 618 CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name="t_orthogonal_converted_to_fm") 619 CALL cp_fm_struct_release(fm_struct_tmp) 620 !CALL copy_dbcsr_to_fm(almo_scf_env%matrix_t(ispin), mo_fm_copy) 621 622 CALL init_mo_set(mos(ispin)%mo_set, fm_ref=mo_fm_copy, name='fm_mo') 623 624 CALL cp_fm_release(mo_fm_copy) 625 626 END DO 627 628 CALL timestop(handle) 629 630 END SUBROUTINE construct_qs_mos 631 632! ************************************************************************************************** 633!> \brief return density matrix to the qs_env 634!> \param qs_env ... 635!> \param matrix_p ... 636!> \param mat_distr_aos ... 637!> \par History 638!> 2011.05 created [Rustam Z Khaliullin] 639!> \author Rustam Z Khaliullin 640! ************************************************************************************************** 641 SUBROUTINE almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos) 642 TYPE(qs_environment_type), POINTER :: qs_env 643 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p 644 INTEGER, INTENT(IN) :: mat_distr_aos 645 646 CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_env', & 647 routineP = moduleN//':'//routineN 648 649 INTEGER :: handle, ispin, nspins 650 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao 651 TYPE(qs_rho_type), POINTER :: rho 652 653 CALL timeset(routineN, handle) 654 655 NULLIFY (rho, rho_ao) 656 nspins = SIZE(matrix_p) 657 CALL get_qs_env(qs_env, rho=rho) 658 CALL qs_rho_get(rho, rho_ao=rho_ao) 659 660 ! set the new density matrix 661 DO ispin = 1, nspins 662 CALL matrix_almo_to_qs(matrix_p(ispin), & 663 rho_ao(ispin)%matrix, & 664 mat_distr_aos) 665 END DO 666 CALL qs_rho_update_rho(rho, qs_env=qs_env) 667 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) 668 669 CALL timestop(handle) 670 671 END SUBROUTINE almo_dm_to_qs_env 672 673! ************************************************************************************************** 674!> \brief uses the ALMO density matrix 675!> to compute KS matrix (inside QS environment) and the new energy 676!> \param qs_env ... 677!> \param matrix_p ... 678!> \param energy_total ... 679!> \param mat_distr_aos ... 680!> \param smear ... 681!> \param kTS_sum ... 682!> \par History 683!> 2011.05 created [Rustam Z Khaliullin] 684!> 2018.09 smearing support [Ruben Staub] 685!> \author Rustam Z Khaliullin 686! ************************************************************************************************** 687 SUBROUTINE almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, smear, kTS_sum) 688 TYPE(qs_environment_type), POINTER :: qs_env 689 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p 690 REAL(KIND=dp) :: energy_total 691 INTEGER, INTENT(IN) :: mat_distr_aos 692 LOGICAL, INTENT(IN), OPTIONAL :: smear 693 REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum 694 695 CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_ks', & 696 routineP = moduleN//':'//routineN 697 698 INTEGER :: handle 699 LOGICAL :: smearing 700 REAL(KIND=dp) :: entropic_term 701 TYPE(qs_energy_type), POINTER :: energy 702 703 CALL timeset(routineN, handle) 704 705 IF (PRESENT(smear)) THEN 706 smearing = smear 707 ELSE 708 smearing = .FALSE. 709 ENDIF 710 711 IF (PRESENT(kTS_sum)) THEN 712 entropic_term = kTS_sum 713 ELSE 714 entropic_term = 0.0_dp 715 ENDIF 716 717 NULLIFY (energy) 718 CALL get_qs_env(qs_env, energy=energy) 719 CALL almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos) 720 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., & 721 print_active=.TRUE.) 722 723 !! Add electronic entropy contribution if smearing is requested 724 !! Previous QS entropy is replaced by the sum of the entropy for each spin 725 IF (smearing) THEN 726 energy%total = energy%total - energy%kTS + entropic_term 727 END IF 728 729 energy_total = energy%total 730 731 CALL timestop(handle) 732 733 END SUBROUTINE almo_dm_to_qs_ks 734 735! ************************************************************************************************** 736!> \brief uses the ALMO density matrix 737!> to compute ALMO KS matrix and the new energy 738!> \param qs_env ... 739!> \param matrix_p ... 740!> \param matrix_ks ... 741!> \param energy_total ... 742!> \param eps_filter ... 743!> \param mat_distr_aos ... 744!> \param smear ... 745!> \param kTS_sum ... 746!> \par History 747!> 2011.05 created [Rustam Z Khaliullin] 748!> 2018.09 smearing support [Ruben Staub] 749!> \author Rustam Z Khaliullin 750! ************************************************************************************************** 751 SUBROUTINE almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, & 752 mat_distr_aos, smear, kTS_sum) 753 754 TYPE(qs_environment_type), POINTER :: qs_env 755 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks 756 REAL(KIND=dp) :: energy_total, eps_filter 757 INTEGER, INTENT(IN) :: mat_distr_aos 758 LOGICAL, INTENT(IN), OPTIONAL :: smear 759 REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum 760 761 CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_almo_ks', & 762 routineP = moduleN//':'//routineN 763 764 INTEGER :: handle, ispin, nspins 765 LOGICAL :: smearing 766 REAL(KIND=dp) :: entropic_term 767 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks 768 769 CALL timeset(routineN, handle) 770 771 IF (PRESENT(smear)) THEN 772 smearing = smear 773 ELSE 774 smearing = .FALSE. 775 ENDIF 776 777 IF (PRESENT(kTS_sum)) THEN 778 entropic_term = kTS_sum 779 ELSE 780 entropic_term = 0.0_dp 781 ENDIF 782 783 ! update KS matrix in the QS env 784 CALL almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, & 785 smear=smearing, & 786 kTS_sum=entropic_term) 787 788 nspins = SIZE(matrix_ks) 789 790 ! get KS matrix from the QS env and convert to the ALMO format 791 CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks) 792 DO ispin = 1, nspins 793 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, & 794 matrix_ks(ispin), & 795 mat_distr_aos, .FALSE.) 796 CALL dbcsr_filter(matrix_ks(ispin), eps_filter) 797 ENDDO 798 799 CALL timestop(handle) 800 801 END SUBROUTINE almo_dm_to_almo_ks 802 803! ************************************************************************************************** 804!> \brief update qs_env total energy 805!> \param qs_env ... 806!> \param energy ... 807!> \param energy_singles_corr ... 808!> \par History 809!> 2013.03 created [Rustam Z Khaliullin] 810!> \author Rustam Z Khaliullin 811! ************************************************************************************************** 812 SUBROUTINE almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr) 813 814 TYPE(qs_environment_type), POINTER :: qs_env 815 REAL(KIND=dp), INTENT(IN), OPTIONAL :: energy, energy_singles_corr 816 817 CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_update_ks_energy', & 818 routineP = moduleN//':'//routineN 819 820 TYPE(qs_energy_type), POINTER :: qs_energy 821 822 CALL get_qs_env(qs_env, energy=qs_energy) 823 824 IF (PRESENT(energy_singles_corr)) THEN 825 qs_energy%singles_corr = energy_singles_corr 826 ELSE 827 qs_energy%singles_corr = 0.0_dp 828 ENDIF 829 830 IF (PRESENT(energy)) THEN 831 qs_energy%total = energy 832 ENDIF 833 834 qs_energy%total = qs_energy%total + qs_energy%singles_corr 835 836 END SUBROUTINE almo_scf_update_ks_energy 837 838! ************************************************************************************************** 839!> \brief Creates the matrix that imposes absolute locality on MOs 840!> \param qs_env ... 841!> \param almo_scf_env ... 842!> \par History 843!> 2011.11 created [Rustam Z. Khaliullin] 844!> \author Rustam Z. Khaliullin 845! ************************************************************************************************** 846 SUBROUTINE almo_scf_construct_quencher(qs_env, almo_scf_env) 847 848 TYPE(qs_environment_type), POINTER :: qs_env 849 TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env 850 851 CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_construct_quencher', & 852 routineP = moduleN//':'//routineN 853 854 CHARACTER :: sym 855 INTEGER :: col, contact_atom_1, contact_atom_2, domain_col, domain_map_local_entries, & 856 domain_row, global_entries, global_list_length, grid1, GroupID, handle, hold, iatom, & 857 iatom2, iblock_col, iblock_row, idomain, idomain2, ientry, igrid, ineig, ineighbor, & 858 iNode, inode2, ipair, ispin, jatom, jatom2, jdomain2, local_list_length, & 859 max_domain_neighbors, max_neig, mynode, nblkcols_tot, nblkrows_tot, nblks, ndomains, & 860 neig_temp, nnode2, nNodes, row, unit_nr 861 INTEGER, ALLOCATABLE, DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, & 862 domain_map_global, domain_map_local, first_atom_of_molecule, global_list, & 863 last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu 864 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: domain_grid, domain_neighbor_list, & 865 domain_neighbor_list_excessive 866 LOGICAL :: already_listed, block_active, & 867 delayed_increment, found, & 868 max_neig_fails, tr 869 REAL(KIND=dp) :: contact1_radius, contact2_radius, & 870 distance, distance_squared, overlap, & 871 r0, r1, s0, s1, trial_distance_squared 872 REAL(KIND=dp), DIMENSION(3) :: rab 873 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_new_block 874 TYPE(cell_type), POINTER :: cell 875 TYPE(cp_logger_type), POINTER :: logger 876 TYPE(dbcsr_distribution_type) :: dist 877 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 878 TYPE(dbcsr_type) :: matrix_s_sym 879 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 880 TYPE(neighbor_list_iterator_p_type), & 881 DIMENSION(:), POINTER :: nl_iterator, nl_iterator2 882 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 883 POINTER :: sab_almo 884 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 885 886 CALL timeset(routineN, handle) 887 888 ! get a useful output_unit 889 logger => cp_get_default_logger() 890 IF (logger%para_env%ionode) THEN 891 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 892 ELSE 893 unit_nr = -1 894 ENDIF 895 896 ndomains = almo_scf_env%ndomains 897 898 CALL get_qs_env(qs_env=qs_env, & 899 particle_set=particle_set, & 900 molecule_set=molecule_set, & 901 cell=cell, & 902 matrix_s=matrix_s, & 903 sab_almo=sab_almo) 904 905 ! if we are dealing with molecules get info about them 906 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. & 907 almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 908 ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules)) 909 ALLOCATE (last_atom_of_molecule(almo_scf_env%nmolecules)) 910 CALL get_molecule_set_info(molecule_set, & 911 mol_to_first_atom=first_atom_of_molecule, & 912 mol_to_last_atom=last_atom_of_molecule) 913 ENDIF 914 915 ! create a symmetrized copy of the ao overlap 916 CALL dbcsr_create(matrix_s_sym, & 917 template=almo_scf_env%matrix_s(1), & 918 matrix_type=dbcsr_type_no_symmetry) 919 CALL dbcsr_get_info(almo_scf_env%matrix_s(1), & 920 matrix_type=sym) 921 IF (sym .EQ. dbcsr_type_no_symmetry) THEN 922 CALL dbcsr_copy(matrix_s_sym, almo_scf_env%matrix_s(1)) 923 ELSE 924 CALL dbcsr_desymmetrize(almo_scf_env%matrix_s(1), & 925 matrix_s_sym) 926 ENDIF 927 928 ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins)) 929 ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins)) 930 931 !DO ispin=1,almo_scf_env%nspins 932 ispin = 1 933 934 ! create the sparsity template for the occupied orbitals 935 CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t(ispin), & 936 matrix_qs=matrix_s(1)%matrix, & 937 almo_scf_env=almo_scf_env, & 938 name_new="T_QUENCHER", & 939 size_keys=(/almo_mat_dim_aobasis, almo_mat_dim_occ/), & 940 symmetry_new=dbcsr_type_no_symmetry, & 941 spin_key=ispin, & 942 init_domains=.FALSE.) 943 944 ! initialize distance quencher 945 CALL dbcsr_work_create(almo_scf_env%quench_t(ispin), & 946 work_mutable=.TRUE.) 947 948 nblkrows_tot = dbcsr_nblkrows_total(almo_scf_env%quench_t(ispin)) 949 nblkcols_tot = dbcsr_nblkcols_total(almo_scf_env%quench_t(ispin)) 950 951 CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist) 952 CALL dbcsr_distribution_get(dist, numnodes=nNodes, group=GroupID, mynode=mynode) 953 954 ! create global atom neighbor list from the local lists 955 ! first, calculate number of local pairs 956 local_list_length = 0 957 CALL neighbor_list_iterator_create(nl_iterator, sab_almo) 958 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 959 ! nnode - total number of neighbors for iatom 960 ! inode - current neighbor count 961 CALL get_iterator_info(nl_iterator, & 962 iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2) 963 !WRITE(*,*) "GET INFO: ",iatom2, jatom2, inode2, nnode2 964 IF (inode2 == 1) THEN 965 local_list_length = local_list_length + nnode2 966 END IF 967 END DO 968 CALL neighbor_list_iterator_release(nl_iterator) 969 970 ! second, extract the local list to an array 971 ALLOCATE (local_list(2*local_list_length)) 972 local_list(:) = 0 973 local_list_length = 0 974 CALL neighbor_list_iterator_create(nl_iterator2, sab_almo) 975 DO WHILE (neighbor_list_iterate(nl_iterator2) == 0) 976 CALL get_iterator_info(nl_iterator2, & 977 iatom=iatom2, jatom=jatom2) 978 local_list(2*local_list_length + 1) = iatom2 979 local_list(2*local_list_length + 2) = jatom2 980 local_list_length = local_list_length + 1 981 ENDDO ! end loop over pairs of atoms 982 CALL neighbor_list_iterator_release(nl_iterator2) 983 984 ! third, communicate local length to the other nodes 985 ALLOCATE (list_length_cpu(nNodes), list_offset_cpu(nNodes)) 986 CALL mp_allgather(2*local_list_length, list_length_cpu, GroupID) 987 988 ! fourth, create a global list 989 list_offset_cpu(1) = 0 990 DO iNode = 2, nNodes 991 list_offset_cpu(iNode) = list_offset_cpu(iNode - 1) + & 992 list_length_cpu(iNode - 1) 993 ENDDO 994 global_list_length = list_offset_cpu(nNodes) + list_length_cpu(nNodes) 995 996 ! fifth, communicate all list data 997 ALLOCATE (global_list(global_list_length)) 998 CALL mp_allgather(local_list, global_list, & 999 list_length_cpu, list_offset_cpu, GroupID) 1000 DEALLOCATE (list_length_cpu, list_offset_cpu) 1001 DEALLOCATE (local_list) 1002 1003 ! calculate maximum number of atoms surrounding the domain 1004 ALLOCATE (current_number_neighbors(almo_scf_env%ndomains)) 1005 current_number_neighbors(:) = 0 1006 global_list_length = global_list_length/2 1007 DO ipair = 1, global_list_length 1008 iatom2 = global_list(2*(ipair - 1) + 1) 1009 jatom2 = global_list(2*(ipair - 1) + 2) 1010 idomain2 = almo_scf_env%domain_index_of_atom(iatom2) 1011 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2) 1012 ! add to the list 1013 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1 1014 ! add j,i with i,j 1015 IF (idomain2 .NE. jdomain2) THEN 1016 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1 1017 ENDIF 1018 ENDDO 1019 max_domain_neighbors = MAXVAL(current_number_neighbors) 1020 1021 ! use the global atom neighbor list to create a global domain neighbor list 1022 ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors)) 1023 current_number_neighbors(:) = 1 1024 DO ipair = 1, ndomains 1025 domain_neighbor_list_excessive(ipair, 1) = ipair 1026 ENDDO 1027 DO ipair = 1, global_list_length 1028 iatom2 = global_list(2*(ipair - 1) + 1) 1029 jatom2 = global_list(2*(ipair - 1) + 2) 1030 idomain2 = almo_scf_env%domain_index_of_atom(iatom2) 1031 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2) 1032 already_listed = .FALSE. 1033 DO ineighbor = 1, current_number_neighbors(idomain2) 1034 IF (domain_neighbor_list_excessive(idomain2, ineighbor) .EQ. jdomain2) THEN 1035 already_listed = .TRUE. 1036 EXIT 1037 ENDIF 1038 ENDDO 1039 IF (.NOT. already_listed) THEN 1040 ! add to the list 1041 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1 1042 domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2 1043 ! add j,i with i,j 1044 IF (idomain2 .NE. jdomain2) THEN 1045 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1 1046 domain_neighbor_list_excessive(jdomain2, current_number_neighbors(jdomain2)) = idomain2 1047 ENDIF 1048 ENDIF 1049 ENDDO ! end loop over pairs of atoms 1050 DEALLOCATE (global_list) 1051 1052 max_domain_neighbors = MAXVAL(current_number_neighbors) 1053 ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors)) 1054 domain_neighbor_list(:, :) = 0 1055 domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors) 1056 DEALLOCATE (domain_neighbor_list_excessive) 1057 1058 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains)) 1059 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2)) 1060 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0 1061 almo_scf_env%domain_map(ispin)%index1(:) = 0 1062 domain_map_local_entries = 0 1063 1064 ! RZK-warning intermediate [0,1] quencher values are ill-defined 1065 ! for molecules (not continuous and conceptually inadequate) 1066 1067 ! O(N) loop over domain pairs 1068 DO row = 1, nblkrows_tot 1069 DO col = 1, current_number_neighbors(row) 1070 tr = .FALSE. 1071 iblock_row = row 1072 iblock_col = domain_neighbor_list(row, col) 1073 CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin), & 1074 iblock_row, iblock_col, hold) 1075 1076 IF (hold .EQ. mynode) THEN 1077 1078 ! Translate indices of distribution blocks to indices of domain blocks 1079 ! Rows are AOs 1080 domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row) 1081 ! Columns are electrons (i.e. MOs) 1082 domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col) 1083 1084 SELECT CASE (almo_scf_env%constraint_type) 1085 CASE (almo_constraint_block_diagonal) 1086 1087 block_active = .FALSE. 1088 ! type of electron groups 1089 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN 1090 1091 ! type of ao domains 1092 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1093 1094 ! ao domains are molecular / electron groups are molecular 1095 IF (domain_row == domain_col) THEN 1096 block_active = .TRUE. 1097 ENDIF 1098 1099 ELSE ! ao domains are atomic 1100 1101 ! ao domains are atomic / electron groups are molecular 1102 CPABORT("Illegal: atomic domains and molecular groups") 1103 1104 ENDIF 1105 1106 ELSE ! electron groups are atomic 1107 1108 ! type of ao domains 1109 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1110 1111 ! ao domains are molecular / electron groups are atomic 1112 CPABORT("Illegal: molecular domains and atomic groups") 1113 1114 ELSE 1115 1116 ! ao domains are atomic / electron groups are atomic 1117 IF (domain_row == domain_col) THEN 1118 block_active = .TRUE. 1119 ENDIF 1120 1121 ENDIF 1122 1123 ENDIF ! end type of electron groups 1124 1125 IF (block_active) THEN 1126 1127 NULLIFY (p_new_block) 1128 CALL dbcsr_reserve_block2d(almo_scf_env%quench_t(ispin), & 1129 iblock_row, iblock_col, p_new_block) 1130 CPASSERT(ASSOCIATED(p_new_block)) 1131 p_new_block(:, :) = 1.0_dp 1132 1133 IF (domain_map_local_entries .GE. max_domain_neighbors*almo_scf_env%ndomains) THEN 1134 CPABORT("weird... max_domain_neighbors is exceeded") 1135 ENDIF 1136 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row 1137 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col 1138 domain_map_local_entries = domain_map_local_entries + 1 1139 1140 ENDIF 1141 1142 CASE (almo_constraint_ao_overlap) 1143 1144 ! type of electron groups 1145 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN 1146 1147 ! type of ao domains 1148 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1149 1150 ! ao domains are molecular / electron groups are molecular 1151 1152 ! compute the maximum overlap between the atoms of the two molecules 1153 CALL dbcsr_get_block_p(matrix_s_sym, & 1154 iblock_row, iblock_col, p_new_block, found) 1155 IF (found) THEN 1156 ! CPErrorMessage(cp_failure_level,routineP,"S block not found") 1157 ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure) 1158 overlap = MAXVAL(ABS(p_new_block)) 1159 ELSE 1160 overlap = 0.0_dp 1161 ENDIF 1162 1163 ELSE ! ao domains are atomic 1164 1165 ! ao domains are atomic / electron groups are molecular 1166 ! overlap_between_atom_and_molecule(atom=domain_row,molecule=domain_col) 1167 CPABORT("atomic domains and molecular groups - NYI") 1168 1169 ENDIF 1170 1171 ELSE ! electron groups are atomic 1172 1173 ! type of ao domains 1174 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1175 1176 ! ao domains are molecular / electron groups are atomic 1177 ! overlap_between_atom_and_molecule(atom=domain_col,molecule=domain_row) 1178 CPABORT("molecular domains and atomic groups - NYI") 1179 1180 ELSE 1181 1182 ! ao domains are atomic / electron groups are atomic 1183 ! compute max overlap between atoms: domain_row and domain_col 1184 CALL dbcsr_get_block_p(matrix_s_sym, & 1185 iblock_row, iblock_col, p_new_block, found) 1186 IF (found) THEN 1187 ! CPErrorMessage(cp_failure_level,routineP,"S block not found") 1188 ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure) 1189 overlap = MAXVAL(ABS(p_new_block)) 1190 ELSE 1191 overlap = 0.0_dp 1192 ENDIF 1193 1194 ENDIF 1195 1196 ENDIF ! end type of electron groups 1197 1198 s0 = -LOG10(ABS(almo_scf_env%quencher_s0)) 1199 s1 = -LOG10(ABS(almo_scf_env%quencher_s1)) 1200 IF (overlap .EQ. 0.0_dp) THEN 1201 overlap = -LOG10(ABS(almo_scf_env%eps_filter)) + 100.0_dp 1202 ELSE 1203 overlap = -LOG10(overlap) 1204 ENDIF 1205 IF (s0 .LT. 0.0_dp) THEN 1206 CPABORT("S0 is less than zero") 1207 ENDIF 1208 IF (s1 .LE. 0.0_dp) THEN 1209 CPABORT("S1 is less than or equal to zero") 1210 ENDIF 1211 IF (s0 .GE. s1) THEN 1212 CPABORT("S0 is greater than or equal to S1") 1213 ENDIF 1214 1215 ! Fill in non-zero blocks if AOs are close to the electron center 1216 IF (overlap .LT. s1) THEN 1217 NULLIFY (p_new_block) 1218 CALL dbcsr_reserve_block2d(almo_scf_env%quench_t(ispin), & 1219 iblock_row, iblock_col, p_new_block) 1220 CPASSERT(ASSOCIATED(p_new_block)) 1221 IF (overlap .LE. s0) THEN 1222 p_new_block(:, :) = 1.0_dp 1223 !WRITE(*,'(A15,2I7,3F8.3,E11.3)') "INTRA-BLOCKS: ",& 1224 ! iblock_col, iblock_row, s0, s1, overlap, p_new_block(1,1) 1225 ELSE 1226 p_new_block(:, :) = 1.0_dp/(1.0_dp + EXP(-(s0 - s1)/(s0 - overlap) - (s0 - s1)/(overlap - s1))) 1227 !WRITE(*,'(A15,2I7,3F8.3,E11.3)') "INTER-BLOCKS: ",& 1228 ! iblock_col, iblock_row, s0, s1, overlap, p_new_block(1,1) 1229 ENDIF 1230 1231 IF (ABS(p_new_block(1, 1)) .GT. ABS(almo_scf_env%eps_filter)) THEN 1232 IF (domain_map_local_entries .GE. max_domain_neighbors*almo_scf_env%ndomains) THEN 1233 CPABORT("weird... max_domain_neighbors is exceeded") 1234 ENDIF 1235 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row 1236 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col 1237 domain_map_local_entries = domain_map_local_entries + 1 1238 ENDIF 1239 1240 ENDIF 1241 1242 CASE (almo_constraint_distance) 1243 1244 ! type of electron groups 1245 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN 1246 1247 ! type of ao domains 1248 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1249 1250 ! ao domains are molecular / electron groups are molecular 1251 1252 ! compute distance between molecules: domain_row and domain_col 1253 ! distance between molecules is defined as the smallest 1254 ! distance among all atom pairs 1255 IF (domain_row == domain_col) THEN 1256 distance = 0.0_dp 1257 contact_atom_1 = first_atom_of_molecule(domain_row) 1258 contact_atom_2 = first_atom_of_molecule(domain_col) 1259 ELSE 1260 distance_squared = 1.0E+100_dp 1261 contact_atom_1 = -1 1262 contact_atom_2 = -1 1263 DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row) 1264 DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col) 1265 rab(:) = pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell) 1266 trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 1267 IF (trial_distance_squared .LT. distance_squared) THEN 1268 distance_squared = trial_distance_squared 1269 contact_atom_1 = iatom 1270 contact_atom_2 = jatom 1271 ENDIF 1272 ENDDO ! jatom 1273 ENDDO ! iatom 1274 CPASSERT(contact_atom_1 .GT. 0) 1275 distance = SQRT(distance_squared) 1276 ENDIF 1277 1278 ELSE ! ao domains are atomic 1279 1280 ! ao domains are atomic / electron groups are molecular 1281 !distance_between_atom_and_molecule(atom=domain_row,molecule=domain_col) 1282 CPABORT("atomic domains and molecular groups - NYI") 1283 1284 ENDIF 1285 1286 ELSE ! electron groups are atomic 1287 1288 ! type of ao domains 1289 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1290 1291 ! ao domains are molecular / electron groups are atomic 1292 !distance_between_atom_and_molecule(atom=domain_col,molecule=domain_row) 1293 CPABORT("molecular domains and atomic groups - NYI") 1294 1295 ELSE 1296 1297 ! ao domains are atomic / electron groups are atomic 1298 ! compute distance between atoms: domain_row and domain_col 1299 rab(:) = pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell) 1300 distance = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 1301 contact_atom_1 = domain_row 1302 contact_atom_2 = domain_col 1303 1304 ENDIF 1305 1306 ENDIF ! end type of electron groups 1307 1308 ! get atomic radii to compute distance cutoff threshold 1309 IF (almo_scf_env%quencher_radius_type == do_bondparm_covalent) THEN 1310 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, & 1311 rcov=contact1_radius) 1312 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, & 1313 rcov=contact2_radius) 1314 ELSE IF (almo_scf_env%quencher_radius_type == do_bondparm_vdw) THEN 1315 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, & 1316 rvdw=contact1_radius) 1317 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, & 1318 rvdw=contact2_radius) 1319 ELSE 1320 CPABORT("Illegal quencher_radius_type") 1321 END IF 1322 contact1_radius = cp_unit_to_cp2k(contact1_radius, "angstrom") 1323 contact2_radius = cp_unit_to_cp2k(contact2_radius, "angstrom") 1324 1325 !RZK-warning the procedure is faulty for molecules: 1326 ! the closest contacts should be found using 1327 ! the element specific radii 1328 1329 ! compute inner and outer cutoff radii 1330 r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius) 1331 !+almo_scf_env%quencher_r0_shift 1332 r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius) 1333 !+almo_scf_env%quencher_r1_shift 1334 1335 IF (r0 .LT. 0.0_dp) THEN 1336 CPABORT("R0 is less than zero") 1337 ENDIF 1338 IF (r1 .LE. 0.0_dp) THEN 1339 CPABORT("R1 is less than or equal to zero") 1340 ENDIF 1341 IF (r0 .GT. r1) THEN 1342 CPABORT("R0 is greater than or equal to R1") 1343 ENDIF 1344 1345 ! Fill in non-zero blocks if AOs are close to the electron center 1346 IF (distance .LT. r1) THEN 1347 NULLIFY (p_new_block) 1348 CALL dbcsr_reserve_block2d(almo_scf_env%quench_t(ispin), & 1349 iblock_row, iblock_col, p_new_block) 1350 CPASSERT(ASSOCIATED(p_new_block)) 1351 IF (distance .LE. r0) THEN 1352 p_new_block(:, :) = 1.0_dp 1353 !WRITE(*,'(A15,2I7,5F8.3,E11.3)') "INTRA-BLOCKS: ",& 1354 ! iblock_col, iblock_row, contact1_radius,& 1355 ! contact2_radius, r0, r1, distance, p_new_block(1,1) 1356 ELSE 1357 ! remove the intermediate values from the quencher temporarily 1358 CPABORT("") 1359 p_new_block(:, :) = 1.0_dp/(1.0_dp + EXP((r1 - r0)/(r0 - distance) + (r1 - r0)/(r1 - distance))) 1360 !WRITE(*,'(A15,2I7,5F8.3,E11.3)') "INTER-BLOCKS: ",& 1361 ! iblock_col, iblock_row, contact1_radius,& 1362 ! contact2_radius, r0, r1, distance, p_new_block(1,1) 1363 ENDIF 1364 1365 IF (ABS(p_new_block(1, 1)) .GT. ABS(almo_scf_env%eps_filter)) THEN 1366 IF (domain_map_local_entries .GE. max_domain_neighbors*almo_scf_env%ndomains) THEN 1367 CPABORT("weird... max_domain_neighbors is exceeded") 1368 ENDIF 1369 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row 1370 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col 1371 domain_map_local_entries = domain_map_local_entries + 1 1372 ENDIF 1373 1374 ENDIF 1375 1376 CASE DEFAULT 1377 CPABORT("Illegal constraint type") 1378 END SELECT 1379 1380 ENDIF ! mynode 1381 1382 ENDDO 1383 ENDDO ! end O(N) loop over pairs 1384 1385 DEALLOCATE (domain_neighbor_list) 1386 DEALLOCATE (current_number_neighbors) 1387 1388 CALL dbcsr_finalize(almo_scf_env%quench_t(ispin)) 1389 !CALL dbcsr_scale(almo_scf_env%quench_t(ispin),& 1390 ! almo_scf_env%envelope_amplitude) 1391 CALL dbcsr_filter(almo_scf_env%quench_t(ispin), & 1392 almo_scf_env%eps_filter) 1393 1394 ! check that both domain_map and quench_t have the same number of entries 1395 nblks = dbcsr_get_num_blocks(almo_scf_env%quench_t(ispin)) 1396 IF (nblks .NE. domain_map_local_entries) THEN 1397 CPABORT("number of blocks is wrong") 1398 ENDIF 1399 1400 ! first, communicate map sizes on the other nodes 1401 ALLOCATE (domain_entries_cpu(nNodes), offset_for_cpu(nNodes)) 1402 CALL mp_allgather(2*domain_map_local_entries, domain_entries_cpu, GroupID) 1403 1404 ! second, create 1405 offset_for_cpu(1) = 0 1406 DO iNode = 2, nNodes 1407 offset_for_cpu(iNode) = offset_for_cpu(iNode - 1) + & 1408 domain_entries_cpu(iNode - 1) 1409 ENDDO 1410 global_entries = offset_for_cpu(nNodes) + domain_entries_cpu(nNodes) 1411 1412 ! communicate all entries 1413 ALLOCATE (domain_map_global(global_entries)) 1414 ALLOCATE (domain_map_local(2*domain_map_local_entries)) 1415 DO ientry = 1, domain_map_local_entries 1416 domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1) 1417 domain_map_local(2*ientry) = almo_scf_env%domain_map(ispin)%pairs(ientry, 2) 1418 ENDDO 1419 CALL mp_allgather(domain_map_local, domain_map_global, & 1420 domain_entries_cpu, offset_for_cpu, GroupID) 1421 DEALLOCATE (domain_entries_cpu, offset_for_cpu) 1422 DEALLOCATE (domain_map_local) 1423 1424 DEALLOCATE (almo_scf_env%domain_map(ispin)%index1) 1425 DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs) 1426 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains)) 1427 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2)) 1428 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0 1429 almo_scf_env%domain_map(ispin)%index1(:) = 0 1430 1431 ! unpack the received data into a local variable 1432 ! since we do not know the maximum global number of neighbors 1433 ! try one. if fails increase the maximum number and try again 1434 ! until it succeeds 1435 max_neig = max_domain_neighbors 1436 max_neig_fails = .TRUE. 1437 max_neig_loop: DO WHILE (max_neig_fails) 1438 ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig)) 1439 domain_grid(:, :) = 0 1440 ! init the number of collected neighbors 1441 domain_grid(:, 0) = 1 1442 ! loop over the records 1443 global_entries = global_entries/2 1444 DO ientry = 1, global_entries 1445 ! get the center 1446 grid1 = domain_map_global(2*ientry) 1447 ! get the neighbor 1448 ineig = domain_map_global(2*(ientry - 1) + 1) 1449 ! check boundaries 1450 IF (domain_grid(grid1, 0) .GT. max_neig) THEN 1451 ! this neighbor will overstep the boundaries 1452 ! stop the trial and increase the max number of neighbors 1453 DEALLOCATE (domain_grid) 1454 max_neig = max_neig*2 1455 CYCLE max_neig_loop 1456 ENDIF 1457 ! for the current center loop over the collected neighbors 1458 ! to insert the current record in a numerical order 1459 delayed_increment = .FALSE. 1460 DO igrid = 1, domain_grid(grid1, 0) 1461 ! compare the current neighbor with that already in the 'book' 1462 IF (ineig .LT. domain_grid(grid1, igrid)) THEN 1463 ! if this one is smaller then insert it here and pick up the one 1464 ! from the book to continue inserting 1465 neig_temp = ineig 1466 ineig = domain_grid(grid1, igrid) 1467 domain_grid(grid1, igrid) = neig_temp 1468 ELSE 1469 IF (domain_grid(grid1, igrid) .EQ. 0) THEN 1470 ! got the empty slot now - insert the record 1471 domain_grid(grid1, igrid) = ineig 1472 ! increase the record counter but do it outside the loop 1473 delayed_increment = .TRUE. 1474 ENDIF 1475 ENDIF 1476 ENDDO 1477 IF (delayed_increment) THEN 1478 domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1 1479 ELSE 1480 ! should not be here - all records must be inserted 1481 CPABORT("all records must be inserted") 1482 ENDIF 1483 ENDDO 1484 max_neig_fails = .FALSE. 1485 ENDDO max_neig_loop 1486 DEALLOCATE (domain_map_global) 1487 1488 ientry = 1 1489 DO idomain = 1, almo_scf_env%ndomains 1490 DO ineig = 1, domain_grid(idomain, 0) - 1 1491 almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig) 1492 almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain 1493 ientry = ientry + 1 1494 ENDDO 1495 almo_scf_env%domain_map(ispin)%index1(idomain) = ientry 1496 ENDDO 1497 DEALLOCATE (domain_grid) 1498 1499 !ENDDO ! ispin 1500 IF (almo_scf_env%nspins .EQ. 2) THEN 1501 CALL dbcsr_copy(almo_scf_env%quench_t(2), & 1502 almo_scf_env%quench_t(1)) 1503 almo_scf_env%domain_map(2)%pairs(:, :) = & 1504 almo_scf_env%domain_map(1)%pairs(:, :) 1505 almo_scf_env%domain_map(2)%index1(:) = & 1506 almo_scf_env%domain_map(1)%index1(:) 1507 ENDIF 1508 1509 CALL dbcsr_release(matrix_s_sym) 1510 1511 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. & 1512 almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN 1513 DEALLOCATE (first_atom_of_molecule) 1514 DEALLOCATE (last_atom_of_molecule) 1515 ENDIF 1516 1517 !mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(& 1518 ! dbcsr_distribution(almo_scf_env%quench_t(ispin)))) 1519 !nblkrows_tot = dbcsr_nblkrows_total(almo_scf_env%quench_t(ispin)) 1520 !nblkcols_tot = dbcsr_nblkcols_total(almo_scf_env%quench_t(ispin)) 1521 !DO row = 1, nblkrows_tot 1522 ! DO col = 1, nblkcols_tot 1523 ! tr = .FALSE. 1524 ! iblock_row = row 1525 ! iblock_col = col 1526 ! CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin),& 1527 ! iblock_row, iblock_col, tr, hold) 1528 ! CALL dbcsr_get_block_p(almo_scf_env%quench_t(ispin),& 1529 ! row, col, p_new_block, found) 1530 ! write(*,*) "RST_NOTE:", mynode, row, col, hold, found 1531 ! ENDDO 1532 !ENDDO 1533 1534 CALL timestop(handle) 1535 1536 END SUBROUTINE almo_scf_construct_quencher 1537 1538! ***************************************************************************** 1539!> \brief Compute matrix W (energy-weighted density matrix) that is needed 1540!> for the evaluation of forces 1541!> \param matrix_w ... 1542!> \param almo_scf_env ... 1543!> \par History 1544!> 2015.03 created [Rustam Z. Khaliullin] 1545!> \author Rustam Z. Khaliullin 1546! ************************************************************************************************** 1547 SUBROUTINE calculate_w_matrix_almo(matrix_w, almo_scf_env) 1548 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w 1549 TYPE(almo_scf_env_type) :: almo_scf_env 1550 1551 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_almo', & 1552 routineP = moduleN//':'//routineN 1553 1554 INTEGER :: handle, ispin 1555 REAL(KIND=dp) :: scaling 1556 TYPE(dbcsr_type) :: tmp_nn1, tmp_no1, tmp_oo1, tmp_oo2 1557 1558 CALL timeset(routineN, handle) 1559 1560 IF (almo_scf_env%nspins == 1) THEN 1561 scaling = 2.0_dp 1562 ELSE 1563 scaling = 1.0_dp 1564 ENDIF 1565 1566 DO ispin = 1, almo_scf_env%nspins 1567 1568 CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), & 1569 matrix_type=dbcsr_type_no_symmetry) 1570 CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), & 1571 matrix_type=dbcsr_type_no_symmetry) 1572 CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), & 1573 matrix_type=dbcsr_type_no_symmetry) 1574 CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), & 1575 matrix_type=dbcsr_type_no_symmetry) 1576 1577 CALL dbcsr_copy(tmp_nn1, almo_scf_env%matrix_ks(ispin)) 1578 ! 1. TMP_NO1=F.T 1579 CALL dbcsr_multiply("N", "N", scaling, tmp_nn1, almo_scf_env%matrix_t(ispin), & 1580 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter) 1581 ! 2. TMP_OO1=T^(tr).TMP_NO1=T^(tr).(FT) 1582 CALL dbcsr_multiply("T", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_no1, & 1583 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter) 1584 ! 3. TMP_OO2=TMP_OO1.siginv=(T^(tr)FT).siginv 1585 CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_oo1, almo_scf_env%matrix_sigma_inv(ispin), & 1586 0.0_dp, tmp_oo2, filter_eps=almo_scf_env%eps_filter) 1587 ! 4. TMP_OO1=siginv.TMP_OO2=siginv.(T^(tr)FTsiginv) 1588 CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_sigma_inv(ispin), tmp_oo2, & 1589 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter) 1590 ! 5. TMP_NO1=T.TMP_OO1.=T.(siginvT^(tr)FTsiginv) 1591 CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_oo1, & 1592 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter) 1593 ! 6. TMP_NN1=TMP_NO1.T^(tr)=(TsiginvT^(tr)FTsiginv).T^(tr)=RFR 1594 CALL dbcsr_multiply("N", "T", 1.0_dp, tmp_no1, almo_scf_env%matrix_t(ispin), & 1595 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter) 1596 CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos) 1597 1598 CALL dbcsr_release(tmp_nn1) 1599 CALL dbcsr_release(tmp_no1) 1600 CALL dbcsr_release(tmp_oo1) 1601 CALL dbcsr_release(tmp_oo2) 1602 1603 ENDDO 1604 1605 CALL timestop(handle) 1606 1607 END SUBROUTINE calculate_w_matrix_almo 1608 1609END MODULE almo_scf_qs 1610 1611