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