1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Module that collects all Coulomb parts of the fock matrix construction 8!> 9!> \author Teodoro Laino (05.2009) [tlaino] - split and module reorganization 10!> \par History 11!> Teodoro Laino (04.2008) [tlaino] - University of Zurich : d-orbitals 12!> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Speed-up 13!> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Periodic SE 14!> Teodoro Laino (05.2009) [tlaino] - Stress Tensor 15!> 16! ************************************************************************************************** 17MODULE se_fock_matrix_coulomb 18 USE atomic_kind_types, ONLY: atomic_kind_type,& 19 get_atomic_kind_set 20 USE atprop_types, ONLY: atprop_type 21 USE cell_types, ONLY: cell_type 22 USE cp_control_types, ONLY: dft_control_type,& 23 semi_empirical_control_type 24 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& 25 dbcsr_deallocate_matrix_set 26 USE cp_log_handling, ONLY: cp_get_default_logger,& 27 cp_logger_type 28 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 29 cp_print_key_unit_nr 30 USE cp_para_types, ONLY: cp_para_env_type 31 USE dbcsr_api, ONLY: dbcsr_add,& 32 dbcsr_distribute,& 33 dbcsr_get_block_diag,& 34 dbcsr_get_block_p,& 35 dbcsr_p_type,& 36 dbcsr_replicate_all,& 37 dbcsr_set,& 38 dbcsr_sum_replicated 39 USE distribution_1d_types, ONLY: distribution_1d_type 40 USE ewald_environment_types, ONLY: ewald_env_get,& 41 ewald_environment_type 42 USE ewald_pw_types, ONLY: ewald_pw_get,& 43 ewald_pw_type 44 USE ewalds_multipole, ONLY: ewald_multipole_evaluate 45 USE fist_neighbor_list_control, ONLY: list_control 46 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type 47 USE input_constants, ONLY: & 48 do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, & 49 do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1, do_se_IS_slater 50 USE input_section_types, ONLY: section_vals_get_subs_vals,& 51 section_vals_type 52 USE kinds, ONLY: dp 53 USE message_passing, ONLY: mp_sum 54 USE multipole_types, ONLY: do_multipole_charge,& 55 do_multipole_dipole,& 56 do_multipole_none,& 57 do_multipole_quadrupole 58 USE particle_types, ONLY: particle_type 59 USE pw_poisson_types, ONLY: do_ewald_ewald 60 USE qs_energy_types, ONLY: qs_energy_type 61 USE qs_environment_types, ONLY: get_qs_env,& 62 qs_environment_type 63 USE qs_force_types, ONLY: qs_force_type 64 USE qs_kind_types, ONLY: get_qs_kind,& 65 qs_kind_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 se_fock_matrix_integrals, ONLY: & 73 dfock2C, dfock2C_r3, dfock2_1el, dfock2_1el_r3, fock2C, fock2C_ew, fock2C_r3, fock2_1el, & 74 fock2_1el_ew, fock2_1el_r3, se_coulomb_ij_interaction 75 USE semi_empirical_int_arrays, ONLY: rij_threshold,& 76 se_orbital_pointer 77 USE semi_empirical_integrals, ONLY: corecore_el,& 78 dcorecore_el 79 USE semi_empirical_mpole_methods, ONLY: quadrupole_sph_to_cart 80 USE semi_empirical_mpole_types, ONLY: nddo_mpole_type,& 81 semi_empirical_mpole_type 82 USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type 83 USE semi_empirical_types, ONLY: get_se_param,& 84 se_int_control_type,& 85 se_taper_type,& 86 semi_empirical_p_type,& 87 semi_empirical_type,& 88 setup_se_int_control_type 89 USE semi_empirical_utils, ONLY: finalize_se_taper,& 90 get_se_type,& 91 initialize_se_taper 92 USE virial_methods, ONLY: virial_pair_force 93 USE virial_types, ONLY: virial_type 94#include "./base/base_uses.f90" 95 96 IMPLICIT NONE 97 PRIVATE 98 99 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_coulomb' 100 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 101 102 PUBLIC :: build_fock_matrix_coulomb, build_fock_matrix_coulomb_lr, & 103 build_fock_matrix_coul_lr_r3 104 105CONTAINS 106 107! ************************************************************************************************** 108!> \brief Construction of the Coulomb part of the Fock matrix 109!> \param qs_env ... 110!> \param ks_matrix ... 111!> \param matrix_p ... 112!> \param energy ... 113!> \param calculate_forces ... 114!> \param store_int_env ... 115!> \author JGH 116! ************************************************************************************************** 117 SUBROUTINE build_fock_matrix_coulomb(qs_env, ks_matrix, matrix_p, energy, calculate_forces, & 118 store_int_env) 119 120 TYPE(qs_environment_type), POINTER :: qs_env 121 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p 122 TYPE(qs_energy_type), POINTER :: energy 123 LOGICAL, INTENT(in) :: calculate_forces 124 TYPE(semi_empirical_si_type), POINTER :: store_int_env 125 126 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coulomb', & 127 routineP = moduleN//':'//routineN 128 129 INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, ispin, itype, jatom, jkind, natom, & 130 natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nspins 131 INTEGER, ALLOCATABLE, DIMENSION(:) :: se_natorb 132 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 133 LOGICAL :: anag, atener, check, defined, found, & 134 switch, use_virial 135 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 136 REAL(KIND=dp) :: delta, dr1, ecore2, ecores 137 REAL(KIND=dp), DIMENSION(2) :: ecab 138 REAL(KIND=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b 139 REAL(KIND=dp), DIMENSION(3) :: force_ab, rij 140 REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b 141 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, ksb_block_a, & 142 ksb_block_b, pa_block_a, pa_block_b, & 143 pb_block_a, pb_block_b 144 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 145 TYPE(atprop_type), POINTER :: atprop 146 TYPE(cell_type), POINTER :: cell 147 TYPE(cp_para_env_type), POINTER :: para_env 148 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_ks, diagmat_p 149 TYPE(dft_control_type), POINTER :: dft_control 150 TYPE(ewald_environment_type), POINTER :: ewald_env 151 TYPE(ewald_pw_type), POINTER :: ewald_pw 152 TYPE(neighbor_list_iterator_p_type), & 153 DIMENSION(:), POINTER :: nl_iterator 154 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 155 POINTER :: sab_se 156 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 157 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 158 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 159 TYPE(se_int_control_type) :: se_int_control 160 TYPE(se_taper_type), POINTER :: se_taper 161 TYPE(semi_empirical_control_type), POINTER :: se_control 162 TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list 163 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 164 TYPE(virial_type), POINTER :: virial 165 166 CALL timeset(routineN, handle) 167 168 NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, & 169 se_control, se_taper, virial, atprop) 170 171 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, & 172 para_env=para_env, sab_se=sab_se, atomic_kind_set=atomic_kind_set, atprop=atprop, & 173 qs_kind_set=qs_kind_set, particle_set=particle_set, virial=virial) 174 175 ! Parameters 176 CALL initialize_se_taper(se_taper, coulomb=.TRUE.) 177 se_control => dft_control%qs_control%se_control 178 anag = se_control%analytical_gradients 179 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 180 atener = atprop%energy 181 182 CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, & 183 do_ewald_gks=se_control%do_ewald_gks, integral_screening=se_control%integral_screening, & 184 shortrange=(se_control%do_ewald .OR. se_control%do_ewald_gks), & 185 max_multipole=se_control%max_multipole, pc_coulomb_int=.FALSE.) 186 IF (se_control%do_ewald_gks) THEN 187 CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw) 188 CALL ewald_env_get(ewald_env, alpha=se_int_control%ewald_gks%alpha) 189 CALL ewald_pw_get(ewald_pw, pw_big_pool=se_int_control%ewald_gks%pw_pool, & 190 dg=se_int_control%ewald_gks%dg) 191 END IF 192 193 nspins = dft_control%nspins 194 CPASSERT(ASSOCIATED(matrix_p)) 195 CPASSERT(SIZE(ks_matrix) > 0) 196 197 nkind = SIZE(atomic_kind_set) 198 IF (calculate_forces) THEN 199 CALL get_qs_env(qs_env=qs_env, force=force) 200 natom = SIZE(particle_set) 201 delta = se_control%delta 202 ! Allocate atom index for kind 203 ALLOCATE (atom_of_kind(natom)) 204 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 205 END IF 206 207 CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins) 208 CALL dbcsr_allocate_matrix_set(diagmat_p, nspins) 209 210 DO ispin = 1, nspins 211 ! Allocate diagonal block matrices 212 ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr 213 CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix) 214 CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix) 215 CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp) 216 CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix) 217 CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix) 218 END DO 219 220 ecore2 = 0.0_dp 221 itype = get_se_type(dft_control%qs_control%method_id) 222 ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind)) 223 DO ikind = 1, nkind 224 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 225 se_kind_list(ikind)%se_param => se_kind_a 226 CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a) 227 se_defined(ikind) = (defined .AND. natorb_a >= 1) 228 se_natorb(ikind) = natorb_a 229 END DO 230 231 CALL neighbor_list_iterator_create(nl_iterator, sab_se) 232 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 233 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij) 234 IF (.NOT. se_defined(ikind)) CYCLE 235 IF (.NOT. se_defined(jkind)) CYCLE 236 se_kind_a => se_kind_list(ikind)%se_param 237 se_kind_b => se_kind_list(jkind)%se_param 238 natorb_a = se_natorb(ikind) 239 natorb_b = se_natorb(jkind) 240 natorb_a2 = natorb_a**2 241 natorb_b2 = natorb_b**2 242 243 IF (inode == 1) THEN 244 CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, & 245 row=iatom, col=iatom, BLOCK=pa_block_a, found=found) 246 CPASSERT(ASSOCIATED(pa_block_a)) 247 check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a) 248 CPASSERT(check) 249 CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, & 250 row=iatom, col=iatom, BLOCK=ksa_block_a, found=found) 251 CPASSERT(ASSOCIATED(ksa_block_a)) 252 p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a 253 pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/)) 254 IF (nspins == 2) THEN 255 CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, & 256 row=iatom, col=iatom, BLOCK=pa_block_b, found=found) 257 CPASSERT(ASSOCIATED(pa_block_b)) 258 check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a) 259 CPASSERT(check) 260 CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, & 261 row=iatom, col=iatom, BLOCK=ksa_block_b, found=found) 262 CPASSERT(ASSOCIATED(ksa_block_b)) 263 p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b 264 pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/)) 265 END IF 266 267 ENDIF 268 269 dr1 = DOT_PRODUCT(rij, rij) 270 IF (dr1 > rij_threshold) THEN 271 ! Determine the order of the atoms, and in case switch them.. 272 IF (iatom <= jatom) THEN 273 switch = .FALSE. 274 ELSE 275 switch = .TRUE. 276 END IF 277 ! Retrieve blocks for KS and P 278 CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, & 279 row=jatom, col=jatom, BLOCK=pb_block_a, found=found) 280 CPASSERT(ASSOCIATED(pb_block_a)) 281 check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b) 282 CPASSERT(check) 283 CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, & 284 row=jatom, col=jatom, BLOCK=ksb_block_a, found=found) 285 CPASSERT(ASSOCIATED(ksb_block_a)) 286 p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a 287 pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/)) 288 ! Handle more than one configuration 289 IF (nspins == 2) THEN 290 CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, & 291 row=jatom, col=jatom, BLOCK=pb_block_b, found=found) 292 CPASSERT(ASSOCIATED(pb_block_b)) 293 check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b) 294 CPASSERT(check) 295 CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, & 296 row=jatom, col=jatom, BLOCK=ksb_block_b, found=found) 297 CPASSERT(ASSOCIATED(ksb_block_b)) 298 check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2)) 299 CPASSERT(check) 300 p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b 301 pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/)) 302 END IF 303 304 SELECT CASE (dft_control%qs_control%method_id) 305 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, & 306 do_method_rm1, do_method_mndod, do_method_pnnl) 307 308 ! Two-centers One-electron terms 309 IF (nspins == 1) THEN 310 ecab = 0._dp 311 CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, & 312 pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, & 313 se_taper=se_taper, store_int_env=store_int_env) 314 ecore2 = ecore2 + ecab(1) + ecab(2) 315 ELSE IF (nspins == 2) THEN 316 ecab = 0._dp 317 CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, & 318 pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag, & 319 se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env) 320 CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_b, ksb_block_b, & 321 pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, & 322 se_taper=se_taper, store_int_env=store_int_env) 323 ecore2 = ecore2 + ecab(1) + ecab(2) 324 END IF 325 IF (atener) THEN 326 atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1) 327 atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2) 328 END IF 329 ! Coulomb Terms 330 IF (nspins == 1) THEN 331 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, & 332 ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 333 store_int_env=store_int_env) 334 ELSE IF (nspins == 2) THEN 335 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, & 336 ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 337 store_int_env=store_int_env) 338 339 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, & 340 ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 341 store_int_env=store_int_env) 342 END IF 343 344 IF (calculate_forces) THEN 345 atom_a = atom_of_kind(iatom) 346 atom_b = atom_of_kind(jatom) 347 348 ! Derivatives of the Two-centre One-electron terms 349 force_ab = 0.0_dp 350 IF (nspins == 1) THEN 351 CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_a, pb_a, itype=itype, anag=anag, & 352 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, & 353 delta=delta) 354 ELSE IF (nspins == 2) THEN 355 CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_block_a, pb_block_a, itype=itype, anag=anag, & 356 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 357 CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_b, pb_b, itype=itype, anag=anag, & 358 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 359 END IF 360 IF (use_virial) THEN 361 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij) 362 END IF 363 364 ! Sum up force components 365 force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1) 366 force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1) 367 368 force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2) 369 force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2) 370 371 force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3) 372 force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3) 373 374 ! Derivatives of the Coulomb Terms 375 force_ab = 0._dp 376 IF (nspins == 1) THEN 377 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, & 378 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 379 ELSE IF (nspins == 2) THEN 380 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, & 381 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 382 383 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, & 384 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 385 END IF 386 IF (switch) THEN 387 force_ab(1) = -force_ab(1) 388 force_ab(2) = -force_ab(2) 389 force_ab(3) = -force_ab(3) 390 END IF 391 IF (use_virial) THEN 392 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij) 393 END IF 394 ! Sum up force components 395 force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1) 396 force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1) 397 398 force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2) 399 force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2) 400 401 force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3) 402 force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3) 403 END IF 404 CASE DEFAULT 405 CPABORT("") 406 END SELECT 407 ELSE 408 IF (se_int_control%do_ewald_gks) THEN 409 CPASSERT(iatom == jatom) 410 ! Two-centers One-electron terms 411 ecores = 0._dp 412 IF (nspins == 1) THEN 413 CALL fock2_1el_ew(se_kind_a, rij, ksa_block_a, pa_a, & 414 ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, & 415 se_taper=se_taper, store_int_env=store_int_env) 416 ELSE IF (nspins == 2) THEN 417 CALL fock2_1el_ew(se_kind_a, rij, ksa_block_a, pa_block_a, & 418 ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, & 419 se_taper=se_taper, store_int_env=store_int_env) 420 CALL fock2_1el_ew(se_kind_a, rij, ksa_block_b, pa_b, & 421 ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, & 422 se_taper=se_taper, store_int_env=store_int_env) 423 END IF 424 ecore2 = ecore2 + ecores 425 IF (atener) THEN 426 atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecores 427 END IF 428 ! Coulomb Terms 429 IF (nspins == 1) THEN 430 CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a, & 431 factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 432 store_int_env=store_int_env) 433 ELSE IF (nspins == 2) THEN 434 CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a, & 435 factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 436 store_int_env=store_int_env) 437 CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_b, & 438 factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 439 store_int_env=store_int_env) 440 END IF 441 END IF 442 END IF 443 END DO 444 CALL neighbor_list_iterator_release(nl_iterator) 445 446 DEALLOCATE (se_kind_list, se_defined, se_natorb) 447 448 DO ispin = 1, nspins 449 CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix) 450 CALL dbcsr_distribute(diagmat_ks(ispin)%matrix) 451 CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, & 452 1.0_dp, 1.0_dp) 453 END DO 454 CALL dbcsr_deallocate_matrix_set(diagmat_p) 455 CALL dbcsr_deallocate_matrix_set(diagmat_ks) 456 457 IF (calculate_forces) THEN 458 DEALLOCATE (atom_of_kind) 459 END IF 460 461 ! Two-centers one-electron terms 462 CALL mp_sum(ecore2, para_env%group) 463 energy%hartree = ecore2 - energy%core 464! WRITE ( *, * ) 'IN SE_F_COUL', ecore2, energy%core 465 466 CALL finalize_se_taper(se_taper) 467 468 CALL timestop(handle) 469 END SUBROUTINE build_fock_matrix_coulomb 470 471! ************************************************************************************************** 472!> \brief Long-Range part for SE Coulomb interactions 473!> \param qs_env ... 474!> \param ks_matrix ... 475!> \param matrix_p ... 476!> \param energy ... 477!> \param calculate_forces ... 478!> \param store_int_env ... 479!> \date 08.2008 [created] 480!> \author Teodoro Laino [tlaino] - University of Zurich 481! ************************************************************************************************** 482 SUBROUTINE build_fock_matrix_coulomb_lr(qs_env, ks_matrix, matrix_p, energy, & 483 calculate_forces, store_int_env) 484 485 TYPE(qs_environment_type), POINTER :: qs_env 486 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p 487 TYPE(qs_energy_type), POINTER :: energy 488 LOGICAL, INTENT(in) :: calculate_forces 489 TYPE(semi_empirical_si_type), POINTER :: store_int_env 490 491 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coulomb_lr', & 492 routineP = moduleN//':'//routineN 493 494 INTEGER :: atom_a, ewald_type, forces_g_size, handle, iatom, ikind, ilist, indi, indj, & 495 ispin, itype, iw, jint, natoms, natorb_a, nkind, nlocal_particles, node, nparticle_local, & 496 nspins, size_1c_int 497 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 498 LOGICAL :: anag, atener, defined, found, use_virial 499 LOGICAL, DIMENSION(3) :: task 500 REAL(KIND=dp) :: e_neut, e_self, energy_glob, & 501 energy_local, enuc, fac, tmp 502 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: forces_g, forces_r 503 REAL(KIND=dp), DIMENSION(3) :: force_a 504 REAL(KIND=dp), DIMENSION(3, 3) :: pv_glob, pv_local, qcart 505 REAL(KIND=dp), DIMENSION(5) :: qsph 506 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, pa_block_a 507 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 508 TYPE(atprop_type), POINTER :: atprop 509 TYPE(cell_type), POINTER :: cell 510 TYPE(cp_logger_type), POINTER :: logger 511 TYPE(cp_para_env_type), POINTER :: para_env 512 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_ks, diagmat_p 513 TYPE(dft_control_type), POINTER :: dft_control 514 TYPE(distribution_1d_type), POINTER :: local_particles 515 TYPE(ewald_environment_type), POINTER :: ewald_env 516 TYPE(ewald_pw_type), POINTER :: ewald_pw 517 TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env 518 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 519 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 520 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 521 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 522 TYPE(section_vals_type), POINTER :: se_section 523 TYPE(semi_empirical_control_type), POINTER :: se_control 524 TYPE(semi_empirical_mpole_type), POINTER :: mpole 525 TYPE(semi_empirical_type), POINTER :: se_kind_a 526 TYPE(virial_type), POINTER :: virial 527 528 CALL timeset(routineN, handle) 529 530 NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, local_particles, & 531 se_control, ewald_env, ewald_pw, se_nddo_mpole, se_nonbond_env, se_section, mpole, & 532 logger, virial, atprop) 533 534 logger => cp_get_default_logger() 535 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env, & 536 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set, & 537 ewald_env=ewald_env, & 538 local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole, & 539 se_nonbond_env=se_nonbond_env, virial=virial, atprop=atprop) 540 541 nlocal_particles = SUM(local_particles%n_el(:)) 542 natoms = SIZE(particle_set) 543 CALL ewald_env_get(ewald_env, ewald_type=ewald_type) 544 SELECT CASE (ewald_type) 545 CASE (do_ewald_ewald) 546 forces_g_size = nlocal_particles 547 CASE DEFAULT 548 CPABORT("Periodic SE implemented only for standard EWALD sums.") 549 END SELECT 550 551 ! Parameters 552 se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE") 553 se_control => dft_control%qs_control%se_control 554 anag = se_control%analytical_gradients 555 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces 556 atener = atprop%energy 557 558 nspins = dft_control%nspins 559 CPASSERT(ASSOCIATED(matrix_p)) 560 CPASSERT(SIZE(ks_matrix) > 0) 561 562 CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins) 563 CALL dbcsr_allocate_matrix_set(diagmat_p, nspins) 564 565 nkind = SIZE(atomic_kind_set) 566 DO ispin = 1, nspins 567 ! Allocate diagonal block matrices 568 ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr 569 CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix) 570 CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix) 571 CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp) 572 CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix) 573 CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix) 574 END DO 575 576 ! Check for implemented SE methods 577 SELECT CASE (dft_control%qs_control%method_id) 578 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, & 579 do_method_rm1, do_method_mndod, do_method_pnnl) 580 itype = get_se_type(dft_control%qs_control%method_id) 581 CASE DEFAULT 582 CPABORT("") 583 END SELECT 584 585 ! Zero arrays and possibly build neighbor lists 586 energy_local = 0.0_dp 587 energy_glob = 0.0_dp 588 e_neut = 0.0_dp 589 e_self = 0.0_dp 590 task = .FALSE. 591 SELECT CASE (se_control%max_multipole) 592 CASE (do_multipole_none) 593 ! Do Nothing 594 CASE (do_multipole_charge) 595 task(1) = .TRUE. 596 CASE (do_multipole_dipole) 597 task = .TRUE. 598 task(3) = .FALSE. 599 CASE (do_multipole_quadrupole) 600 task = .TRUE. 601 CASE DEFAULT 602 CPABORT("") 603 END SELECT 604 605 ! Build-up neighbor lists for real-space part of Ewald multipoles 606 CALL list_control(atomic_kind_set, particle_set, local_particles, & 607 cell, se_nonbond_env, para_env, se_section) 608 609 enuc = 0.0_dp 610 energy%core_overlap = 0.0_dp 611 se_nddo_mpole%charge = 0.0_dp 612 se_nddo_mpole%dipole = 0.0_dp 613 se_nddo_mpole%quadrupole = 0.0_dp 614 615 DO ispin = 1, nspins 616 ! Compute the NDDO mpole expansion 617 DO ikind = 1, nkind 618 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 619 CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a) 620 621 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 622 623 nparticle_local = local_particles%n_el(ikind) 624 DO ilist = 1, nparticle_local 625 iatom = local_particles%list(ikind)%array(ilist) 626 CALL dbcsr_get_block_p(matrix=diagmat_p(ispin)%matrix, & 627 row=iatom, col=iatom, BLOCK=pa_block_a, found=found) 628 CPASSERT(ASSOCIATED(pa_block_a)) 629 ! Nuclei 630 IF (task(1) .AND. ispin == 1) se_nddo_mpole%charge(iatom) = se_kind_a%zeff 631 ! Electrons 632 size_1c_int = SIZE(se_kind_a%w_mpole) 633 DO jint = 1, size_1c_int 634 mpole => se_kind_a%w_mpole(jint)%mpole 635 indi = se_orbital_pointer(mpole%indi) 636 indj = se_orbital_pointer(mpole%indj) 637 fac = 1.0_dp 638 IF (indi /= indj) fac = 2.0_dp 639 640 ! Charge 641 IF (mpole%task(1) .AND. task(1)) THEN 642 se_nddo_mpole%charge(iatom) = se_nddo_mpole%charge(iatom) + & 643 fac*pa_block_a(indi, indj)*mpole%c 644 END IF 645 646 ! Dipole 647 IF (mpole%task(2) .AND. task(2)) THEN 648 se_nddo_mpole%dipole(:, iatom) = se_nddo_mpole%dipole(:, iatom) + & 649 fac*pa_block_a(indi, indj)*mpole%d(:) 650 END IF 651 652 ! Quadrupole 653 IF (mpole%task(3) .AND. task(3)) THEN 654 qsph = fac*mpole%qs*pa_block_a(indi, indj) 655 CALL quadrupole_sph_to_cart(qcart, qsph) 656 se_nddo_mpole%quadrupole(:, :, iatom) = se_nddo_mpole%quadrupole(:, :, iatom) + & 657 qcart 658 END IF 659 END DO 660 ! Print some info about charge, dipole and quadrupole (debug purpose only) 661 IF (debug_this_module) THEN 662 WRITE (*, '(I5,F12.6,5X,3F12.6,5X,9F12.6)') iatom, se_nddo_mpole%charge(iatom), & 663 se_nddo_mpole%dipole(:, iatom), se_nddo_mpole%quadrupole(:, :, iatom) 664 END IF 665 END DO 666 END DO 667 END DO 668 CALL mp_sum(se_nddo_mpole%charge, para_env%group) 669 CALL mp_sum(se_nddo_mpole%dipole, para_env%group) 670 CALL mp_sum(se_nddo_mpole%quadrupole, para_env%group) 671 672 ! Initialize for virial 673 IF (use_virial) THEN 674 pv_glob = 0.0_dp 675 pv_local = 0.0_dp 676 END IF 677 678 ! Ewald Multipoles Sum 679 iw = cp_print_key_unit_nr(logger, se_section, "PRINT%EWALD_INFO", extension=".seLog") 680 IF (calculate_forces) THEN 681 CALL get_qs_env(qs_env=qs_env, force=force) 682 683 ! Allocate atom index for kind 684 ALLOCATE (atom_of_kind(natoms)) 685 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 686 687 ! Allocate and zeroing arrays 688 ALLOCATE (forces_g(3, forces_g_size)) 689 ALLOCATE (forces_r(3, natoms)) 690 forces_g = 0.0_dp 691 forces_r = 0.0_dp 692 CALL ewald_multipole_evaluate( & 693 ewald_env, ewald_pw, se_nonbond_env, cell, & 694 particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, & 695 do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=use_virial, do_efield=.TRUE., & 696 charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, & 697 forces_local=forces_g, forces_glob=forces_r, pv_glob=pv_glob, pv_local=pv_local, & 698 efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, iw=iw, & 699 do_debug=.TRUE.) 700 ! Only SR force have to be summed up.. the one in g-space are already fully local.. 701 CALL mp_sum(forces_r, para_env%group) 702 ELSE 703 CALL ewald_multipole_evaluate( & 704 ewald_env, ewald_pw, se_nonbond_env, cell, & 705 particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, & 706 do_correction_bonded=.FALSE., do_forces=.FALSE., do_stress=.FALSE., do_efield=.TRUE., & 707 charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, & 708 efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, & 709 iw=iw, do_debug=.TRUE.) 710 END IF 711 CALL cp_print_key_finished_output(iw, logger, se_section, "PRINT%EWALD_INFO") 712 713 ! Apply correction only when the Integral Scheme is different from Slater 714 IF ((se_control%integral_screening /= do_se_IS_slater) .AND. (.NOT. debug_this_module)) THEN 715 CALL build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, calculate_forces, & 716 store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, virial, & 717 pv_glob) 718 END IF 719 720 ! Virial for the long-range part and correction 721 IF (use_virial) THEN 722 ! Sum up contribution of pv_glob on each thread and keep only one copy of pv_local 723 virial%pv_virial = virial%pv_virial + pv_glob 724 IF (para_env%ionode) THEN 725 virial%pv_virial = virial%pv_virial + pv_local 726 END IF 727 END IF 728 729 ! Debug Statements 730 IF (debug_this_module) THEN 731 CALL mp_sum(energy_glob, para_env%group) 732 WRITE (*, *) "TOTAL ENERGY AFTER EWALD:", energy_local + energy_glob + e_neut + e_self, & 733 energy_local, energy_glob, e_neut, e_self 734 END IF 735 736 ! Modify the KS matrix and possibly compute derivatives 737 node = 0 738 DO ikind = 1, nkind 739 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 740 CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a) 741 742 IF (.NOT. defined .OR. natorb_a < 1) CYCLE 743 744 nparticle_local = local_particles%n_el(ikind) 745 DO ilist = 1, nparticle_local 746 node = node + 1 747 iatom = local_particles%list(ikind)%array(ilist) 748 DO ispin = 1, nspins 749 CALL dbcsr_get_block_p(matrix=diagmat_ks(ispin)%matrix, & 750 row=iatom, col=iatom, BLOCK=ksa_block_a, found=found) 751 CPASSERT(ASSOCIATED(ksa_block_a)) 752 753 ! Modify Hamiltonian Matrix accordingly potential, field and electric field gradient 754 size_1c_int = SIZE(se_kind_a%w_mpole) 755 DO jint = 1, size_1c_int 756 tmp = 0.0_dp 757 mpole => se_kind_a%w_mpole(jint)%mpole 758 indi = se_orbital_pointer(mpole%indi) 759 indj = se_orbital_pointer(mpole%indj) 760 761 ! Charge 762 IF (mpole%task(1) .AND. task(1)) THEN 763 tmp = tmp + mpole%c*se_nddo_mpole%efield0(iatom) 764 END IF 765 766 ! Dipole 767 IF (mpole%task(2) .AND. task(2)) THEN 768 tmp = tmp - DOT_PRODUCT(mpole%d, se_nddo_mpole%efield1(:, iatom)) 769 END IF 770 771 ! Quadrupole 772 IF (mpole%task(3) .AND. task(3)) THEN 773 tmp = tmp - (1.0_dp/3.0_dp)*SUM(mpole%qc*RESHAPE(se_nddo_mpole%efield2(:, iatom), (/3, 3/))) 774 END IF 775 776 ksa_block_a(indi, indj) = ksa_block_a(indi, indj) + tmp 777 ksa_block_a(indj, indi) = ksa_block_a(indi, indj) 778 END DO 779 END DO 780 781 ! Nuclear term and forces 782 IF (task(1)) enuc = enuc + se_kind_a%zeff*se_nddo_mpole%efield0(iatom) 783 IF (atener) THEN 784 atprop%atecoul(iatom) = atprop%atecoul(iatom) + & 785 0.5_dp*se_kind_a%zeff*se_nddo_mpole%efield0(iatom) 786 END IF 787 IF (calculate_forces) THEN 788 atom_a = atom_of_kind(iatom) 789 force_a = forces_r(1:3, iatom) + forces_g(1:3, node) 790 ! Derivatives of the periodic Coulomb Terms 791 force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_a(:) 792 END IF 793 END DO 794 END DO 795 ! Sum nuclear energy contribution 796 CALL mp_sum(enuc, para_env%group) 797 energy%core_overlap = energy%core_overlap + energy%core_overlap0 + 0.5_dp*enuc 798 799 ! Debug Statements 800 IF (debug_this_module) THEN 801 WRITE (*, *) "ENUC: ", enuc*0.5_dp 802 END IF 803 804 DO ispin = 1, nspins 805 CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix) 806 CALL dbcsr_distribute(diagmat_ks(ispin)%matrix) 807 CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, & 808 1.0_dp, 1.0_dp) 809 END DO 810 CALL dbcsr_deallocate_matrix_set(diagmat_p) 811 CALL dbcsr_deallocate_matrix_set(diagmat_ks) 812 813 ! Set the Fock matrix contribution to SCP 814 IF (calculate_forces) THEN 815 DEALLOCATE (atom_of_kind) 816 DEALLOCATE (forces_g) 817 DEALLOCATE (forces_r) 818 END IF 819 820 CALL timestop(handle) 821 END SUBROUTINE build_fock_matrix_coulomb_lr 822 823! ************************************************************************************************** 824!> \brief When doing long-range SE calculation, this module computes the correction 825!> between the mismatch of point-like multipoles and multipoles represented 826!> with charges 827!> \param qs_env ... 828!> \param ks_matrix ... 829!> \param matrix_p ... 830!> \param energy ... 831!> \param calculate_forces ... 832!> \param store_int_env ... 833!> \param se_nddo_mpole ... 834!> \param task ... 835!> \param diagmat_p ... 836!> \param diagmat_ks ... 837!> \param virial ... 838!> \param pv_glob ... 839!> \author Teodoro Laino [tlaino] - 05.2009 840! ************************************************************************************************** 841 SUBROUTINE build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, & 842 calculate_forces, store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, & 843 virial, pv_glob) 844 845 TYPE(qs_environment_type), POINTER :: qs_env 846 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p 847 TYPE(qs_energy_type), POINTER :: energy 848 LOGICAL, INTENT(in) :: calculate_forces 849 TYPE(semi_empirical_si_type), POINTER :: store_int_env 850 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 851 LOGICAL, DIMENSION(3), INTENT(IN) :: task 852 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_p, diagmat_ks 853 TYPE(virial_type), POINTER :: virial 854 REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv_glob 855 856 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coul_lrc', & 857 routineP = moduleN//':'//routineN 858 859 INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, itype, jatom, jkind, natom, & 860 natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nspins, size1, size2 861 INTEGER, ALLOCATABLE, DIMENSION(:) :: se_natorb 862 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 863 LOGICAL :: anag, atener, check, defined, found, & 864 switch, use_virial 865 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 866 REAL(KIND=dp) :: delta, dr1, ecore2, enuc, enuclear, & 867 ptens11, ptens12, ptens13, ptens21, & 868 ptens22, ptens23, ptens31, ptens32, & 869 ptens33 870 REAL(KIND=dp), DIMENSION(2) :: ecab 871 REAL(KIND=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b 872 REAL(KIND=dp), DIMENSION(3) :: force_ab, force_ab0, rij 873 REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b 874 REAL(KIND=dp), DIMENSION(:), POINTER :: efield0 875 REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2, ksa_block_a, ksa_block_b, & 876 ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, pb_block_a, pb_block_b 877 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 878 TYPE(atprop_type), POINTER :: atprop 879 TYPE(cell_type), POINTER :: cell 880 TYPE(cp_para_env_type), POINTER :: para_env 881 TYPE(dft_control_type), POINTER :: dft_control 882 TYPE(neighbor_list_iterator_p_type), & 883 DIMENSION(:), POINTER :: nl_iterator 884 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 885 POINTER :: sab_lrc 886 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 887 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 888 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 889 TYPE(se_int_control_type) :: se_int_control 890 TYPE(se_taper_type), POINTER :: se_taper 891 TYPE(semi_empirical_control_type), POINTER :: se_control 892 TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list 893 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 894 895 CALL timeset(routineN, handle) 896 NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper, & 897 efield0, efield1, efield2, atprop) 898 899 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, & 900 para_env=para_env, sab_lrc=sab_lrc, atomic_kind_set=atomic_kind_set, & 901 qs_kind_set=qs_kind_set, particle_set=particle_set, atprop=atprop) 902 903 ! Parameters 904 CALL initialize_se_taper(se_taper, lr_corr=.TRUE.) 905 se_control => dft_control%qs_control%se_control 906 anag = se_control%analytical_gradients 907 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces 908 atener = atprop%energy 909 910 CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, & 911 do_ewald_gks=.FALSE., integral_screening=se_control%integral_screening, & 912 shortrange=.FALSE., max_multipole=se_control%max_multipole, & 913 pc_coulomb_int=.TRUE.) 914 915 nspins = dft_control%nspins 916 CPASSERT(ASSOCIATED(matrix_p)) 917 CPASSERT(SIZE(ks_matrix) > 0) 918 CPASSERT(ASSOCIATED(diagmat_p)) 919 CPASSERT(ASSOCIATED(diagmat_ks)) 920 921 nkind = SIZE(atomic_kind_set) 922 IF (calculate_forces) THEN 923 CALL get_qs_env(qs_env=qs_env, force=force) 924 natom = SIZE(particle_set) 925 delta = se_control%delta 926 ! Allocate atom index for kind 927 ALLOCATE (atom_of_kind(natom)) 928 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 929 END IF 930 931 ! Allocate arrays for storing partial information on potential, field, field gradient 932 size1 = SIZE(se_nddo_mpole%efield0) 933 ALLOCATE (efield0(size1)) 934 efield0 = 0.0_dp 935 size1 = SIZE(se_nddo_mpole%efield1, 1) 936 size2 = SIZE(se_nddo_mpole%efield1, 2) 937 ALLOCATE (efield1(size1, size2)) 938 efield1 = 0.0_dp 939 size1 = SIZE(se_nddo_mpole%efield2, 1) 940 size2 = SIZE(se_nddo_mpole%efield2, 2) 941 ALLOCATE (efield2(size1, size2)) 942 efield2 = 0.0_dp 943 944 ! Initialize if virial is requested 945 IF (use_virial) THEN 946 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp 947 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp 948 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp 949 END IF 950 951 ! Start of the loop for the correction of the pair interactions 952 ecore2 = 0.0_dp 953 enuclear = 0.0_dp 954 itype = get_se_type(dft_control%qs_control%method_id) 955 956 ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind)) 957 DO ikind = 1, nkind 958 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 959 se_kind_list(ikind)%se_param => se_kind_a 960 CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a) 961 se_defined(ikind) = (defined .AND. natorb_a >= 1) 962 se_natorb(ikind) = natorb_a 963 END DO 964 965 CALL neighbor_list_iterator_create(nl_iterator, sab_lrc) 966 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 967 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij) 968 IF (.NOT. se_defined(ikind)) CYCLE 969 IF (.NOT. se_defined(jkind)) CYCLE 970 se_kind_a => se_kind_list(ikind)%se_param 971 se_kind_b => se_kind_list(jkind)%se_param 972 natorb_a = se_natorb(ikind) 973 natorb_b = se_natorb(jkind) 974 natorb_a2 = natorb_a**2 975 natorb_b2 = natorb_b**2 976 977 IF (inode == 1) THEN 978 CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, & 979 row=iatom, col=iatom, BLOCK=pa_block_a, found=found) 980 CPASSERT(ASSOCIATED(pa_block_a)) 981 check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a) 982 CPASSERT(check) 983 CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, & 984 row=iatom, col=iatom, BLOCK=ksa_block_a, found=found) 985 CPASSERT(ASSOCIATED(ksa_block_a)) 986 p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a 987 pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/)) 988 IF (nspins == 2) THEN 989 CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, & 990 row=iatom, col=iatom, BLOCK=pa_block_b, found=found) 991 CPASSERT(ASSOCIATED(pa_block_b)) 992 check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a) 993 CPASSERT(check) 994 CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, & 995 row=iatom, col=iatom, BLOCK=ksa_block_b, found=found) 996 CPASSERT(ASSOCIATED(ksa_block_b)) 997 p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b 998 pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/)) 999 END IF 1000 1001 ENDIF 1002 1003 dr1 = DOT_PRODUCT(rij, rij) 1004 IF (dr1 > rij_threshold) THEN 1005 ! Determine the order of the atoms, and in case switch them.. 1006 IF (iatom <= jatom) THEN 1007 switch = .FALSE. 1008 ELSE 1009 switch = .TRUE. 1010 END IF 1011 1012 ! Point-like interaction corrections 1013 CALL se_coulomb_ij_interaction(iatom, jatom, task, do_forces=calculate_forces, & 1014 do_efield=.TRUE., do_stress=use_virial, charges=se_nddo_mpole%charge, & 1015 dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, & 1016 force_ab=force_ab0, efield0=efield0, efield1=efield1, efield2=efield2, & 1017 rab2=dr1, rab=rij, ptens11=ptens11, ptens12=ptens12, ptens13=ptens13, & 1018 ptens21=ptens21, ptens22=ptens22, ptens23=ptens23, ptens31=ptens31, & 1019 ptens32=ptens32, ptens33=ptens33) 1020 1021 ! Retrieve blocks for KS and P 1022 CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, & 1023 row=jatom, col=jatom, BLOCK=pb_block_a, found=found) 1024 CPASSERT(ASSOCIATED(pb_block_a)) 1025 check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b) 1026 CPASSERT(check) 1027 CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, & 1028 row=jatom, col=jatom, BLOCK=ksb_block_a, found=found) 1029 CPASSERT(ASSOCIATED(ksb_block_a)) 1030 p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a 1031 pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/)) 1032 ! Handle more than one configuration 1033 IF (nspins == 2) THEN 1034 CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, & 1035 row=jatom, col=jatom, BLOCK=pb_block_b, found=found) 1036 CPASSERT(ASSOCIATED(pb_block_b)) 1037 check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b) 1038 CPASSERT(check) 1039 CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, & 1040 row=jatom, col=jatom, BLOCK=ksb_block_b, found=found) 1041 CPASSERT(ASSOCIATED(ksb_block_b)) 1042 check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2)) 1043 CPASSERT(check) 1044 p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b 1045 pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/)) 1046 END IF 1047 1048 SELECT CASE (dft_control%qs_control%method_id) 1049 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, & 1050 do_method_rm1, do_method_mndod) 1051 ! Evaluate nuclear contribution.. 1052 CALL corecore_el(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, & 1053 se_int_control=se_int_control, se_taper=se_taper) 1054 enuclear = enuclear + enuc 1055 1056 ! Two-centers One-electron terms 1057 IF (nspins == 1) THEN 1058 ecab = 0._dp 1059 CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, & 1060 pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, & 1061 se_taper=se_taper, store_int_env=store_int_env) 1062 ecore2 = ecore2 + ecab(1) + ecab(2) 1063 ELSE IF (nspins == 2) THEN 1064 ecab = 0._dp 1065 CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, & 1066 pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag, & 1067 se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env) 1068 CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_b, ksb_block_b, & 1069 pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, & 1070 se_taper=se_taper, store_int_env=store_int_env) 1071 ecore2 = ecore2 + ecab(1) + ecab(2) 1072 END IF 1073 IF (atener) THEN 1074 atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1) + 0.5_dp*enuc 1075 atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2) + 0.5_dp*enuc 1076 END IF 1077 ! Coulomb Terms 1078 IF (nspins == 1) THEN 1079 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, & 1080 ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 1081 store_int_env=store_int_env) 1082 ELSE IF (nspins == 2) THEN 1083 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, & 1084 ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 1085 store_int_env=store_int_env) 1086 1087 CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, & 1088 ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, & 1089 store_int_env=store_int_env) 1090 END IF 1091 1092 IF (calculate_forces) THEN 1093 atom_a = atom_of_kind(iatom) 1094 atom_b = atom_of_kind(jatom) 1095 1096 ! Evaluate nuclear contribution.. 1097 CALL dcorecore_el(se_kind_a, se_kind_b, rij, denuc=force_ab, itype=itype, delta=delta, & 1098 anag=anag, se_int_control=se_int_control, se_taper=se_taper) 1099 1100 ! Derivatives of the Two-centre One-electron terms 1101 IF (nspins == 1) THEN 1102 CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_a, pb_a, itype=itype, anag=anag, & 1103 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, & 1104 delta=delta) 1105 ELSE IF (nspins == 2) THEN 1106 CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_block_a, pb_block_a, itype=itype, anag=anag, & 1107 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 1108 CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_b, pb_b, itype=itype, anag=anag, & 1109 se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 1110 END IF 1111 IF (use_virial) THEN 1112 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij) 1113 END IF 1114 force_ab = force_ab + force_ab0 1115 1116 ! Sum up force components 1117 force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1) 1118 force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1) 1119 1120 force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2) 1121 force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2) 1122 1123 force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3) 1124 force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3) 1125 1126 ! Derivatives of the Coulomb Terms 1127 force_ab = 0._dp 1128 IF (nspins == 1) THEN 1129 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, & 1130 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 1131 ELSE IF (nspins == 2) THEN 1132 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, & 1133 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 1134 1135 CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, & 1136 anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta) 1137 END IF 1138 IF (switch) THEN 1139 force_ab(1) = -force_ab(1) 1140 force_ab(2) = -force_ab(2) 1141 force_ab(3) = -force_ab(3) 1142 END IF 1143 IF (use_virial) THEN 1144 CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij) 1145 END IF 1146 1147 ! Sum up force components 1148 force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1) 1149 force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1) 1150 1151 force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2) 1152 force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2) 1153 1154 force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3) 1155 force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3) 1156 END IF 1157 CASE DEFAULT 1158 CPABORT("") 1159 END SELECT 1160 END IF 1161 END DO 1162 CALL neighbor_list_iterator_release(nl_iterator) 1163 1164 DEALLOCATE (se_kind_list, se_defined, se_natorb) 1165 1166 ! Sum-up Virial constribution (long-range correction) 1167 IF (use_virial) THEN 1168 pv_glob(1, 1) = pv_glob(1, 1) + ptens11 1169 pv_glob(1, 2) = pv_glob(1, 2) + (ptens12 + ptens21)*0.5_dp 1170 pv_glob(1, 3) = pv_glob(1, 3) + (ptens13 + ptens31)*0.5_dp 1171 pv_glob(2, 1) = pv_glob(1, 2) 1172 pv_glob(2, 2) = pv_glob(2, 2) + ptens22 1173 pv_glob(2, 3) = pv_glob(2, 3) + (ptens23 + ptens32)*0.5_dp 1174 pv_glob(3, 1) = pv_glob(1, 3) 1175 pv_glob(3, 2) = pv_glob(2, 3) 1176 pv_glob(3, 3) = pv_glob(3, 3) + ptens33 1177 END IF 1178 1179 IF (calculate_forces) THEN 1180 DEALLOCATE (atom_of_kind) 1181 END IF 1182 1183 ! collect information on potential, field, field gradient 1184 CALL mp_sum(efield0, para_env%group) 1185 CALL mp_sum(efield1, para_env%group) 1186 CALL mp_sum(efield2, para_env%group) 1187 se_nddo_mpole%efield0 = se_nddo_mpole%efield0 - efield0 1188 se_nddo_mpole%efield1 = se_nddo_mpole%efield1 - efield1 1189 se_nddo_mpole%efield2 = se_nddo_mpole%efield2 - efield2 1190 ! deallocate working arrays 1191 DEALLOCATE (efield0) 1192 DEALLOCATE (efield1) 1193 DEALLOCATE (efield2) 1194 1195 ! Corrections for two-centers one-electron terms + nuclear terms 1196 CALL mp_sum(enuclear, para_env%group) 1197 CALL mp_sum(ecore2, para_env%group) 1198 energy%hartree = energy%hartree + ecore2 1199 energy%core_overlap = enuclear 1200 CALL finalize_se_taper(se_taper) 1201 CALL timestop(handle) 1202 END SUBROUTINE build_fock_matrix_coul_lrc 1203 1204! ************************************************************************************************** 1205!> \brief Construction of the residual part (1/R^3) of the Coulomb long-range 1206!> term of the Fock matrix 1207!> The 1/R^3 correction works in real-space strictly on the zero-cell, 1208!> in order to avoid more parameters to be provided in the input.. 1209!> \param qs_env ... 1210!> \param ks_matrix ... 1211!> \param matrix_p ... 1212!> \param energy ... 1213!> \param calculate_forces ... 1214!> \author Teodoro Laino [tlaino] - 12.2008 1215! ************************************************************************************************** 1216 SUBROUTINE build_fock_matrix_coul_lr_r3(qs_env, ks_matrix, matrix_p, energy, calculate_forces) 1217 TYPE(qs_environment_type), POINTER :: qs_env 1218 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p 1219 TYPE(qs_energy_type), POINTER :: energy 1220 LOGICAL, INTENT(in) :: calculate_forces 1221 1222 CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coul_lr_r3', & 1223 routineP = moduleN//':'//routineN 1224 1225 INTEGER :: atom_a, atom_b, ewald_type, handle, iatom, ikind, inode, ispin, itype, jatom, & 1226 jkind, natoms, natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nlocal_particles, nspins 1227 INTEGER, ALLOCATABLE, DIMENSION(:) :: se_natorb 1228 INTEGER, DIMENSION(:), POINTER :: atom_of_kind 1229 LOGICAL :: anag, atener, check, defined, found, & 1230 switch 1231 LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined 1232 REAL(KIND=dp) :: dr1, ecore2, r2inv, r3inv, rinv 1233 REAL(KIND=dp), DIMENSION(2) :: ecab 1234 REAL(KIND=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b 1235 REAL(KIND=dp), DIMENSION(3) :: dr3inv, force_ab, rij 1236 REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b 1237 REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, ksb_block_a, & 1238 ksb_block_b, pa_block_a, pa_block_b, & 1239 pb_block_a, pb_block_b 1240 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1241 TYPE(atprop_type), POINTER :: atprop 1242 TYPE(cell_type), POINTER :: cell 1243 TYPE(cp_logger_type), POINTER :: logger 1244 TYPE(cp_para_env_type), POINTER :: para_env 1245 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_ks, diagmat_p 1246 TYPE(dft_control_type), POINTER :: dft_control 1247 TYPE(distribution_1d_type), POINTER :: local_particles 1248 TYPE(ewald_environment_type), POINTER :: ewald_env 1249 TYPE(ewald_pw_type), POINTER :: ewald_pw 1250 TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env 1251 TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole 1252 TYPE(neighbor_list_iterator_p_type), & 1253 DIMENSION(:), POINTER :: nl_iterator 1254 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 1255 POINTER :: sab_orb 1256 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1257 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 1258 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 1259 TYPE(section_vals_type), POINTER :: se_section 1260 TYPE(semi_empirical_control_type), POINTER :: se_control 1261 TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list 1262 TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b 1263 1264 CALL timeset(routineN, handle) 1265 1266 CALL get_qs_env(qs_env=qs_env) !sm->dbcsr 1267 1268 NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, & 1269 diagmat_p, local_particles, se_control, ewald_env, ewald_pw, & 1270 se_nddo_mpole, se_nonbond_env, se_section, sab_orb, logger, atprop) 1271 1272 logger => cp_get_default_logger() 1273 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env, & 1274 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set, & 1275 ewald_env=ewald_env, atprop=atprop, & 1276 local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole, & 1277 se_nonbond_env=se_nonbond_env, sab_orb=sab_orb) 1278 1279 nlocal_particles = SUM(local_particles%n_el(:)) 1280 natoms = SIZE(particle_set) 1281 CALL ewald_env_get(ewald_env, ewald_type=ewald_type) 1282 SELECT CASE (ewald_type) 1283 CASE (do_ewald_ewald) 1284 ! Do Nothing 1285 CASE DEFAULT 1286 CPABORT("Periodic SE implemented only for standard EWALD sums.") 1287 END SELECT 1288 1289 ! Parameters 1290 se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE") 1291 se_control => dft_control%qs_control%se_control 1292 anag = se_control%analytical_gradients 1293 atener = atprop%energy 1294 1295 nspins = dft_control%nspins 1296 CPASSERT(ASSOCIATED(matrix_p)) 1297 CPASSERT(SIZE(ks_matrix) > 0) 1298 1299 CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins) 1300 CALL dbcsr_allocate_matrix_set(diagmat_p, nspins) 1301 1302 nkind = SIZE(atomic_kind_set) 1303 DO ispin = 1, nspins 1304 ! Allocate diagonal block matrices 1305 ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr 1306 CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix) 1307 CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix) 1308 CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp) 1309 CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix) 1310 CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix) 1311 END DO 1312 1313 ! Possibly compute forces 1314 IF (calculate_forces) THEN 1315 CALL get_qs_env(qs_env=qs_env, force=force) 1316 ALLOCATE (atom_of_kind(natoms)) 1317 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 1318 END IF 1319 itype = get_se_type(dft_control%qs_control%method_id) 1320 1321 ecore2 = 0.0_dp 1322 ! Real space part of the 1/R^3 sum 1323 ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind)) 1324 DO ikind = 1, nkind 1325 CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a) 1326 se_kind_list(ikind)%se_param => se_kind_a 1327 CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a) 1328 se_defined(ikind) = (defined .AND. natorb_a >= 1) 1329 se_natorb(ikind) = natorb_a 1330 END DO 1331 1332 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 1333 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 1334 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij) 1335 IF (.NOT. se_defined(ikind)) CYCLE 1336 IF (.NOT. se_defined(jkind)) CYCLE 1337 se_kind_a => se_kind_list(ikind)%se_param 1338 se_kind_b => se_kind_list(jkind)%se_param 1339 natorb_a = se_natorb(ikind) 1340 natorb_b = se_natorb(jkind) 1341 natorb_a2 = natorb_a**2 1342 natorb_b2 = natorb_b**2 1343 1344 IF (inode == 1) THEN 1345 CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, & 1346 row=iatom, col=iatom, BLOCK=pa_block_a, found=found) 1347 CPASSERT(ASSOCIATED(pa_block_a)) 1348 check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a) 1349 CPASSERT(check) 1350 CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, & 1351 row=iatom, col=iatom, BLOCK=ksa_block_a, found=found) 1352 CPASSERT(ASSOCIATED(ksa_block_a)) 1353 p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a 1354 pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/)) 1355 IF (nspins == 2) THEN 1356 CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, & 1357 row=iatom, col=iatom, BLOCK=pa_block_b, found=found) 1358 CPASSERT(ASSOCIATED(pa_block_b)) 1359 check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a) 1360 CPASSERT(check) 1361 CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, & 1362 row=iatom, col=iatom, BLOCK=ksa_block_b, found=found) 1363 CPASSERT(ASSOCIATED(ksa_block_b)) 1364 p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b 1365 pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/)) 1366 END IF 1367 END IF 1368 1369 dr1 = DOT_PRODUCT(rij, rij) 1370 IF (dr1 > rij_threshold) THEN 1371 ! Determine the order of the atoms, and in case switch them.. 1372 IF (iatom <= jatom) THEN 1373 switch = .FALSE. 1374 ELSE 1375 switch = .TRUE. 1376 END IF 1377 ! Retrieve blocks for KS and P 1378 CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, & 1379 row=jatom, col=jatom, BLOCK=pb_block_a, found=found) 1380 CPASSERT(ASSOCIATED(pb_block_a)) 1381 check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b) 1382 CPASSERT(check) 1383 CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, & 1384 row=jatom, col=jatom, BLOCK=ksb_block_a, found=found) 1385 CPASSERT(ASSOCIATED(ksb_block_a)) 1386 p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a 1387 pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/)) 1388 ! Handle more than one configuration 1389 IF (nspins == 2) THEN 1390 CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, & 1391 row=jatom, col=jatom, BLOCK=pb_block_b, found=found) 1392 CPASSERT(ASSOCIATED(pb_block_b)) 1393 check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b) 1394 CPASSERT(check) 1395 CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, & 1396 row=jatom, col=jatom, BLOCK=ksb_block_b, found=found) 1397 CPASSERT(ASSOCIATED(ksb_block_b)) 1398 check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2)) 1399 CPASSERT(check) 1400 p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b 1401 pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/)) 1402 END IF 1403 1404 SELECT CASE (dft_control%qs_control%method_id) 1405 CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, & 1406 do_method_rm1, do_method_mndod, do_method_pnnl) 1407 1408 ! Pre-compute some quantities.. 1409 r2inv = 1.0_dp/dr1 1410 rinv = SQRT(r2inv) 1411 r3inv = rinv**3 1412 1413 ! Two-centers One-electron terms 1414 IF (nspins == 1) THEN 1415 ecab = 0._dp 1416 CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_a, pb_a, & 1417 ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, & 1418 e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv) 1419 ecore2 = ecore2 + ecab(1) + ecab(2) 1420 ELSE IF (nspins == 2) THEN 1421 ecab = 0._dp 1422 CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_block_a, & 1423 pb_block_a, ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, & 1424 e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv) 1425 1426 CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_b, ksb_block_b, pa_b, pb_b, & 1427 ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, & 1428 e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv) 1429 ecore2 = ecore2 + ecab(1) + ecab(2) 1430 END IF 1431 IF (atener) THEN 1432 atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1) 1433 atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2) 1434 END IF 1435 ! Coulomb Terms 1436 IF (nspins == 1) THEN 1437 CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, & 1438 ksb_block_a, factor=0.5_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv) 1439 ELSE IF (nspins == 2) THEN 1440 CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, & 1441 ksb_block_a, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv) 1442 1443 CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, & 1444 ksb_block_b, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv) 1445 END IF 1446 1447 ! Compute forces if requested 1448 IF (calculate_forces) THEN 1449 dr3inv = -3.0_dp*rij*r3inv*r2inv 1450 atom_a = atom_of_kind(iatom) 1451 atom_b = atom_of_kind(jatom) 1452 1453 force_ab = 0.0_dp 1454 ! Derivatives of the One-centre One-electron terms 1455 IF (nspins == 1) THEN 1456 CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_a, pb_a, force_ab, & 1457 se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a) 1458 ELSE IF (nspins == 2) THEN 1459 CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_block_a, pb_block_a, force_ab, & 1460 se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a) 1461 1462 CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_b, pb_b, force_ab, & 1463 se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a) 1464 END IF 1465 1466 ! Sum up force components 1467 force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1) 1468 force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1) 1469 1470 force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2) 1471 force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2) 1472 1473 force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3) 1474 force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3) 1475 1476 ! Derivatives of the Coulomb Terms 1477 force_ab = 0.0_dp 1478 IF (nspins == 1) THEN 1479 CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, & 1480 w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab) 1481 ELSE IF (nspins == 2) THEN 1482 CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, & 1483 w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab) 1484 1485 CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, & 1486 w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab) 1487 END IF 1488 1489 ! Sum up force components 1490 force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1) 1491 force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1) 1492 1493 force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2) 1494 force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2) 1495 1496 force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3) 1497 force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3) 1498 END IF 1499 CASE DEFAULT 1500 CPABORT("") 1501 END SELECT 1502 END IF 1503 END DO 1504 CALL neighbor_list_iterator_release(nl_iterator) 1505 1506 DEALLOCATE (se_kind_list, se_defined, se_natorb) 1507 1508 DO ispin = 1, nspins 1509 CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix) 1510 CALL dbcsr_distribute(diagmat_ks(ispin)%matrix) 1511 CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, & 1512 1.0_dp, 1.0_dp) 1513 END DO 1514 CALL dbcsr_deallocate_matrix_set(diagmat_p) 1515 CALL dbcsr_deallocate_matrix_set(diagmat_ks) 1516 1517 IF (calculate_forces) THEN 1518 DEALLOCATE (atom_of_kind) 1519 END IF 1520 1521 ! Two-centers one-electron terms 1522 CALL mp_sum(ecore2, para_env%group) 1523 energy%hartree = energy%hartree + ecore2 1524 1525 CALL timestop(handle) 1526 END SUBROUTINE build_fock_matrix_coul_lr_r3 1527 1528END MODULE se_fock_matrix_coulomb 1529 1530