1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculates the energy contribution and the mo_derivative of 8!> a static electric field (nonperiodic) 9!> \par History 10!> Adjusted from qs_efield_local 11!> \author JGH (10.2019) 12! ************************************************************************************************** 13MODULE ec_efield_local 14 USE ai_moments, ONLY: dipole_force 15 USE atomic_kind_types, ONLY: atomic_kind_type,& 16 get_atomic_kind,& 17 get_atomic_kind_set 18 USE basis_set_types, ONLY: gto_basis_set_p_type,& 19 gto_basis_set_type 20 USE cell_types, ONLY: cell_type,& 21 pbc 22 USE cp_control_types, ONLY: dft_control_type 23 USE cp_para_types, ONLY: cp_para_env_type 24 USE dbcsr_api, ONLY: dbcsr_add,& 25 dbcsr_copy,& 26 dbcsr_get_block_p,& 27 dbcsr_p_type,& 28 dbcsr_set 29 USE ec_env_types, ONLY: energy_correction_type 30 USE kinds, ONLY: dp 31 USE orbital_pointers, ONLY: ncoset 32 USE particle_types, ONLY: particle_type 33 USE qs_energy_types, ONLY: qs_energy_type 34 USE qs_environment_types, ONLY: get_qs_env,& 35 qs_environment_type 36 USE qs_force_types, ONLY: qs_force_type 37 USE qs_kind_types, ONLY: get_qs_kind,& 38 qs_kind_type 39 USE qs_moments, ONLY: build_local_moment_matrix 40 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 41 neighbor_list_iterate,& 42 neighbor_list_iterator_create,& 43 neighbor_list_iterator_p_type,& 44 neighbor_list_iterator_release,& 45 neighbor_list_set_p_type 46 USE qs_period_efield_types, ONLY: efield_berry_type,& 47 init_efield_matrices,& 48 set_efield_matrices 49#include "./base/base_uses.f90" 50 51 IMPLICIT NONE 52 53 PRIVATE 54 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_efield_local' 56 57 ! *** Public subroutines *** 58 59 PUBLIC :: ec_efield_local_operator, ec_efield_integrals 60 61! ************************************************************************************************** 62 63CONTAINS 64 65! ************************************************************************************************** 66 67! ************************************************************************************************** 68!> \brief ... 69!> \param qs_env ... 70!> \param ec_env ... 71!> \param calculate_forces ... 72! ************************************************************************************************** 73 SUBROUTINE ec_efield_local_operator(qs_env, ec_env, calculate_forces) 74 75 TYPE(qs_environment_type), POINTER :: qs_env 76 TYPE(energy_correction_type), POINTER :: ec_env 77 LOGICAL, INTENT(IN) :: calculate_forces 78 79 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_local_operator', & 80 routineP = moduleN//':'//routineN 81 82 INTEGER :: handle 83 REAL(dp), DIMENSION(3) :: rpoint 84 TYPE(dft_control_type), POINTER :: dft_control 85 86 CALL timeset(routineN, handle) 87 88 NULLIFY (dft_control) 89 CALL get_qs_env(qs_env, dft_control=dft_control) 90 91 IF (dft_control%apply_efield) THEN 92 rpoint = 0.0_dp 93 CALL ec_efield_integrals(qs_env, ec_env, rpoint) 94 CALL ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces) 95 END IF 96 97 CALL timestop(handle) 98 99 END SUBROUTINE ec_efield_local_operator 100 101! ************************************************************************************************** 102!> \brief ... 103!> \param qs_env ... 104!> \param ec_env ... 105!> \param rpoint ... 106! ************************************************************************************************** 107 SUBROUTINE ec_efield_integrals(qs_env, ec_env, rpoint) 108 109 TYPE(qs_environment_type), POINTER :: qs_env 110 TYPE(energy_correction_type), POINTER :: ec_env 111 REAL(dp), DIMENSION(3), INTENT(IN) :: rpoint 112 113 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_integrals', & 114 routineP = moduleN//':'//routineN 115 116 INTEGER :: handle, i 117 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s 118 TYPE(dft_control_type), POINTER :: dft_control 119 TYPE(efield_berry_type), POINTER :: efield, efieldref 120 121 CALL timeset(routineN, handle) 122 123 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, efield=efieldref) 124 efield => ec_env%efield 125 CALL init_efield_matrices(efield) 126 matrix_s => ec_env%matrix_s(:, 1) 127 ALLOCATE (dipmat(3)) 128 DO i = 1, 3 129 ALLOCATE (dipmat(i)%matrix) 130 CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'DIP MAT') 131 CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp) 132 END DO 133 CALL build_local_moment_matrix(qs_env, dipmat, 1, rpoint, basis_type="HARRIS") 134 CALL set_efield_matrices(efield=efield, dipmat=dipmat) 135 ec_env%efield => efield 136 137 CALL timestop(handle) 138 139 END SUBROUTINE ec_efield_integrals 140 141! ************************************************************************************************** 142!> \brief ... 143!> \param qs_env ... 144!> \param ec_env ... 145!> \param rpoint ... 146!> \param calculate_forces ... 147! ************************************************************************************************** 148 SUBROUTINE ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces) 149 TYPE(qs_environment_type), POINTER :: qs_env 150 TYPE(energy_correction_type), POINTER :: ec_env 151 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rpoint 152 LOGICAL :: calculate_forces 153 154 CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_mo_derivatives', & 155 routineP = moduleN//':'//routineN 156 157 INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, & 158 jatom, jkind, jset, ldab, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb 159 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 160 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 161 npgfb, nsgfa, nsgfb 162 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 163 LOGICAL :: found, trans 164 REAL(dp) :: charge, dab, fdir 165 REAL(dp), DIMENSION(3) :: ci, fieldpol, ra, rab, rac, rbc, ria 166 REAL(dp), DIMENSION(3, 3) :: forcea, forceb 167 REAL(dp), DIMENSION(:, :), POINTER :: p_block_a, p_block_b, pblock, pmat, work 168 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 169 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 170 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 171 TYPE(cell_type), POINTER :: cell 172 TYPE(cp_para_env_type), POINTER :: para_env 173 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_ks 174 TYPE(dft_control_type), POINTER :: dft_control 175 TYPE(efield_berry_type), POINTER :: efield 176 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 177 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 178 TYPE(neighbor_list_iterator_p_type), & 179 DIMENSION(:), POINTER :: nl_iterator 180 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 181 POINTER :: sab_orb 182 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 183 TYPE(qs_energy_type), POINTER :: energy 184 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 185 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 186 TYPE(qs_kind_type), POINTER :: qs_kind 187 188 CALL timeset(routineN, handle) 189 190 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set) 191 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 192 energy=energy, para_env=para_env, sab_orb=sab_orb) 193 194 efield => ec_env%efield 195 196 fieldpol = dft_control%efield_fields(1)%efield%polarisation* & 197 dft_control%efield_fields(1)%efield%strength 198 199 ! nuclear contribution 200 natom = SIZE(particle_set) 201 IF (calculate_forces) THEN 202 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force) 203 ALLOCATE (atom_of_kind(natom)) 204 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 205 END IF 206 ci = 0.0_dp 207 DO ia = 1, natom 208 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind) 209 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge) 210 ria = particle_set(ia)%r - rpoint 211 ria = pbc(ria, cell) 212 ci(:) = ci(:) + charge*ria(:) 213 IF (calculate_forces) THEN 214 IF (para_env%mepos == 0) THEN 215 iatom = atom_of_kind(ia) 216 DO idir = 1, 3 217 force(ikind)%efield(idir, iatom) = force(ikind)%efield(idir, iatom) - fieldpol(idir)*charge 218 END DO 219 END IF 220 END IF 221 END DO 222 223 IF (ec_env%should_update) THEN 224 ec_env%efield_nuclear = -SUM(ci(:)*fieldpol(:)) 225 226 ! Update KS matrix 227 matrix_ks => ec_env%matrix_h(:, 1) 228 dipmat => efield%dipmat 229 DO ispin = 1, SIZE(matrix_ks) 230 DO idir = 1, 3 231 CALL dbcsr_add(matrix_ks(ispin)%matrix, dipmat(idir)%matrix, & 232 alpha_scalar=1.0_dp, beta_scalar=fieldpol(idir)) 233 END DO 234 END DO 235 END IF 236 237 ! forces from the efield contribution 238 IF (calculate_forces) THEN 239 nkind = SIZE(qs_kind_set) 240 natom = SIZE(particle_set) 241 242 ALLOCATE (basis_set_list(nkind)) 243 DO ikind = 1, nkind 244 qs_kind => qs_kind_set(ikind) 245 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type="HARRIS") 246 IF (ASSOCIATED(basis_set_a)) THEN 247 basis_set_list(ikind)%gto_basis_set => basis_set_a 248 ELSE 249 NULLIFY (basis_set_list(ikind)%gto_basis_set) 250 END IF 251 END DO 252 ! 253 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 254 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 255 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 256 iatom=iatom, jatom=jatom, r=rab) 257 basis_set_a => basis_set_list(ikind)%gto_basis_set 258 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 259 basis_set_b => basis_set_list(jkind)%gto_basis_set 260 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 261 ! basis ikind 262 first_sgfa => basis_set_a%first_sgf 263 la_max => basis_set_a%lmax 264 la_min => basis_set_a%lmin 265 npgfa => basis_set_a%npgf 266 nseta = basis_set_a%nset 267 nsgfa => basis_set_a%nsgf_set 268 rpgfa => basis_set_a%pgf_radius 269 set_radius_a => basis_set_a%set_radius 270 sphi_a => basis_set_a%sphi 271 zeta => basis_set_a%zet 272 ! basis jkind 273 first_sgfb => basis_set_b%first_sgf 274 lb_max => basis_set_b%lmax 275 lb_min => basis_set_b%lmin 276 npgfb => basis_set_b%npgf 277 nsetb = basis_set_b%nset 278 nsgfb => basis_set_b%nsgf_set 279 rpgfb => basis_set_b%pgf_radius 280 set_radius_b => basis_set_b%set_radius 281 sphi_b => basis_set_b%sphi 282 zetb => basis_set_b%zet 283 284 atom_a = atom_of_kind(iatom) 285 atom_b = atom_of_kind(jatom) 286 287 ra(:) = particle_set(iatom)%r(:) - rpoint(:) 288 rac(:) = pbc(ra(:), cell) 289 rbc(:) = rac(:) + rab(:) 290 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 291 292 IF (iatom <= jatom) THEN 293 irow = iatom 294 icol = jatom 295 trans = .FALSE. 296 ELSE 297 irow = jatom 298 icol = iatom 299 trans = .TRUE. 300 END IF 301 302 fdir = 2.0_dp 303 IF (iatom == jatom .AND. dab < 1.e-10_dp) fdir = 1.0_dp 304 305 ! density matrix 306 NULLIFY (p_block_a) 307 CALL dbcsr_get_block_p(ec_env%matrix_p(1, 1)%matrix, irow, icol, p_block_a, found) 308!deb IF (.NOT. found) CYCLE 309 IF (SIZE(ec_env%matrix_p, 1) > 1) THEN 310 NULLIFY (p_block_b) 311 CALL dbcsr_get_block_p(ec_env%matrix_p(2, 1)%matrix, irow, icol, p_block_b, found) 312!deb CPASSERT(found) 313 END IF 314 forcea = 0.0_dp 315 forceb = 0.0_dp 316 317 DO iset = 1, nseta 318 ncoa = npgfa(iset)*ncoset(la_max(iset)) 319 sgfa = first_sgfa(1, iset) 320 DO jset = 1, nsetb 321 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 322 ncob = npgfb(jset)*ncoset(lb_max(jset)) 323 sgfb = first_sgfb(1, jset) 324 ! Calculate the primitive integrals (da|O|b) and (a|O|db) 325 ldab = MAX(ncoa, ncob) 326 ALLOCATE (work(ldab, ldab), pmat(ncoa, ncob)) 327 ! Decontract P matrix block 328 pmat = 0.0_dp 329 DO i = 1, SIZE(ec_env%matrix_p, 1) 330 IF (i == 1) THEN 331 pblock => p_block_a 332 ELSE 333 pblock => p_block_b 334 END IF 335 IF (.NOT. ASSOCIATED(pblock)) CYCLE 336 IF (trans) THEN 337 CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), & 338 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 339 pblock(sgfb, sgfa), SIZE(pblock, 1), & 340 0.0_dp, work(1, 1), ldab) 341 ELSE 342 CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), & 343 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 344 pblock(sgfa, sgfb), SIZE(pblock, 1), & 345 0.0_dp, work(1, 1), ldab) 346 END IF 347 CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), & 348 1.0_dp, work(1, 1), ldab, & 349 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 350 1.0_dp, pmat(1, 1), ncoa) 351 END DO 352 353 CALL dipole_force(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 354 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 355 1, rac, rbc, pmat, forcea, forceb) 356 357 DEALLOCATE (work, pmat) 358 END DO 359 END DO 360 361 DO idir = 1, 3 362 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) & 363 + fdir*fieldpol(idir)*forcea(idir, 1:3) 364 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) & 365 + fdir*fieldpol(idir)*forceb(idir, 1:3) 366 END DO 367 368 END DO 369 CALL neighbor_list_iterator_release(nl_iterator) 370 DEALLOCATE (basis_set_list) 371 DEALLOCATE (atom_of_kind) 372 END IF 373 374 CALL timestop(handle) 375 376 END SUBROUTINE ec_efield_mo_derivatives 377 378END MODULE ec_efield_local 379