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 periodic electric field 9!> \par History 10!> none 11!> \author fschiff (06.2010) 12! ************************************************************************************************** 13MODULE qs_efield_berry 14 USE ai_moments, ONLY: cossin 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 block_p_types, ONLY: block_p_type 21 USE cell_types, ONLY: cell_type,& 22 pbc 23 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm,& 24 cp_cfm_solve 25 USE cp_cfm_types, ONLY: cp_cfm_create,& 26 cp_cfm_p_type,& 27 cp_cfm_release,& 28 cp_cfm_set_all 29 USE cp_control_types, ONLY: dft_control_type 30 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 31 copy_fm_to_dbcsr,& 32 cp_dbcsr_plus_fm_fm_t,& 33 cp_dbcsr_sm_fm_multiply,& 34 dbcsr_deallocate_matrix_set 35 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add 36 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 37 cp_fm_struct_release,& 38 cp_fm_struct_type 39 USE cp_fm_types, ONLY: cp_fm_create,& 40 cp_fm_p_type,& 41 cp_fm_release,& 42 cp_fm_set_all,& 43 cp_fm_type 44 USE cp_gemm_interface, ONLY: cp_gemm 45 USE cp_para_types, ONLY: cp_para_env_type 46 USE dbcsr_api, ONLY: dbcsr_copy,& 47 dbcsr_get_block_p,& 48 dbcsr_p_type,& 49 dbcsr_set,& 50 dbcsr_type 51 USE kinds, ONLY: dp 52 USE mathconstants, ONLY: gaussi,& 53 pi,& 54 twopi,& 55 z_one,& 56 z_zero 57 USE message_passing, ONLY: mp_sum 58 USE orbital_pointers, ONLY: ncoset 59 USE particle_types, ONLY: particle_type 60 USE qs_energy_types, ONLY: qs_energy_type 61 USE qs_environment_types, ONLY: get_qs_env,& 62 qs_environment_type,& 63 set_qs_env 64 USE qs_force_types, ONLY: qs_force_type 65 USE qs_kind_types, ONLY: get_qs_kind,& 66 get_qs_kind_set,& 67 qs_kind_type 68 USE qs_mo_types, ONLY: get_mo_set,& 69 mo_set_p_type 70 USE qs_moments, ONLY: build_berry_moment_matrix 71 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 72 neighbor_list_iterate,& 73 neighbor_list_iterator_create,& 74 neighbor_list_iterator_p_type,& 75 neighbor_list_iterator_release,& 76 neighbor_list_set_p_type 77 USE qs_period_efield_types, ONLY: efield_berry_type,& 78 init_efield_matrices,& 79 set_efield_matrices 80 USE virial_methods, ONLY: virial_pair_force 81 USE virial_types, ONLY: virial_type 82#include "./base/base_uses.f90" 83 84 IMPLICIT NONE 85 86 PRIVATE 87 88 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_efield_berry' 89 90 ! *** Public subroutines *** 91 92 PUBLIC :: qs_efield_berry_phase 93 94! ************************************************************************************************** 95 96CONTAINS 97 98! ************************************************************************************************** 99 100! ************************************************************************************************** 101!> \brief ... 102!> \param qs_env ... 103!> \param just_energy ... 104!> \param calculate_forces ... 105! ************************************************************************************************** 106 SUBROUTINE qs_efield_berry_phase(qs_env, just_energy, calculate_forces) 107 108 TYPE(qs_environment_type), POINTER :: qs_env 109 LOGICAL, INTENT(IN) :: just_energy, calculate_forces 110 111 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_berry_phase', & 112 routineP = moduleN//':'//routineN 113 114 INTEGER :: handle 115 LOGICAL :: s_mstruct_changed 116 TYPE(dft_control_type), POINTER :: dft_control 117 118 CALL timeset(routineN, handle) 119 120 NULLIFY (dft_control) 121 CALL get_qs_env(qs_env, s_mstruct_changed=s_mstruct_changed, & 122 dft_control=dft_control) 123 124 IF (dft_control%apply_period_efield) THEN 125 IF (s_mstruct_changed) CALL qs_efield_integrals(qs_env) 126 IF (dft_control%period_efield%displacement_field) THEN 127 CALL qs_dispfield_derivatives(qs_env, just_energy, calculate_forces) 128 ELSE 129 CALL qs_efield_derivatives(qs_env, just_energy, calculate_forces) 130 END IF 131 END IF 132 133 CALL timestop(handle) 134 135 END SUBROUTINE qs_efield_berry_phase 136 137! ************************************************************************************************** 138!> \brief ... 139!> \param qs_env ... 140! ************************************************************************************************** 141 SUBROUTINE qs_efield_integrals(qs_env) 142 143 TYPE(qs_environment_type), POINTER :: qs_env 144 145 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_integrals', & 146 routineP = moduleN//':'//routineN 147 148 INTEGER :: handle, i 149 REAL(dp), DIMENSION(3) :: kvec 150 TYPE(cell_type), POINTER :: cell 151 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: cosmat, matrix_s, sinmat 152 TYPE(dft_control_type), POINTER :: dft_control 153 TYPE(efield_berry_type), POINTER :: efield 154 155 CALL timeset(routineN, handle) 156 CPASSERT(ASSOCIATED(qs_env)) 157 158 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 159 NULLIFY (matrix_s) 160 CALL get_qs_env(qs_env=qs_env, efield=efield, cell=cell, matrix_s=matrix_s) 161 CALL init_efield_matrices(efield) 162 ALLOCATE (cosmat(3), sinmat(3)) 163 DO i = 1, 3 164 ALLOCATE (cosmat(i)%matrix, sinmat(i)%matrix) 165 166 CALL dbcsr_copy(cosmat(i)%matrix, matrix_s(1)%matrix, 'COS MAT') 167 CALL dbcsr_copy(sinmat(i)%matrix, matrix_s(1)%matrix, 'SIN MAT') 168 CALL dbcsr_set(cosmat(i)%matrix, 0.0_dp) 169 CALL dbcsr_set(sinmat(i)%matrix, 0.0_dp) 170 171 kvec(:) = twopi*cell%h_inv(i, :) 172 CALL build_berry_moment_matrix(qs_env, cosmat(i)%matrix, sinmat(i)%matrix, kvec) 173 END DO 174 CALL set_efield_matrices(efield=efield, cosmat=cosmat, sinmat=sinmat) 175 CALL set_qs_env(qs_env=qs_env, efield=efield) 176 CALL timestop(handle) 177 178 END SUBROUTINE qs_efield_integrals 179 180! ************************************************************************************************** 181!> \brief ... 182!> \param qs_env ... 183!> \param just_energy ... 184!> \param calculate_forces ... 185! ************************************************************************************************** 186 SUBROUTINE qs_efield_derivatives(qs_env, just_energy, calculate_forces) 187 TYPE(qs_environment_type), POINTER :: qs_env 188 LOGICAL, INTENT(IN) :: just_energy, calculate_forces 189 190 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efield_derivatives', & 191 routineP = moduleN//':'//routineN 192 193 COMPLEX(dp) :: zdet, zdeta, zi(3) 194 INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, j, & 195 jatom, jkind, jset, ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, & 196 nseta, nsetb, sgfa, sgfb 197 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 198 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 199 npgfb, nsgfa, nsgfb 200 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 201 LOGICAL :: found, uniform, use_virial 202 REAL(dp) :: charge, ci(3), cqi(3), dab, dd, & 203 ener_field, f0, fab, fieldpol(3), & 204 focc, fpolvec(3), hmat(3, 3), occ, & 205 qi(3), ti(3) 206 REAL(dp), DIMENSION(3) :: forcea, forceb, kvec, ra, rab, rb, ria 207 REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work 208 REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab 209 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 210 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 211 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 212 TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint 213 TYPE(cell_type), POINTER :: cell 214 TYPE(cp_cfm_p_type), DIMENSION(:), POINTER :: eigrmat, inv_mat 215 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_coeff_tmp, mo_derivs_tmp 216 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: inv_work, op_fm_set, opvec 217 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 218 TYPE(cp_fm_type), POINTER :: mo_coeff 219 TYPE(cp_para_env_type), POINTER :: para_env 220 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs 221 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat 222 TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat 223 TYPE(dft_control_type), POINTER :: dft_control 224 TYPE(efield_berry_type), POINTER :: efield 225 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 226 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 227 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 228 TYPE(neighbor_list_iterator_p_type), & 229 DIMENSION(:), POINTER :: nl_iterator 230 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 231 POINTER :: sab_orb 232 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 233 TYPE(qs_energy_type), POINTER :: energy 234 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 235 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 236 TYPE(qs_kind_type), POINTER :: qs_kind 237 TYPE(virial_type), POINTER :: virial 238 239 CALL timeset(routineN, handle) 240 241 NULLIFY (dft_control, cell, particle_set) 242 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, & 243 particle_set=particle_set, virial=virial) 244 NULLIFY (qs_kind_set, efield, para_env, sab_orb) 245 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 246 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb) 247 248 ! calculate stress only if forces requested also 249 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 250 use_virial = use_virial .AND. calculate_forces 251 ! disable stress calculation 252 IF (use_virial) THEN 253 CPABORT("Stress tensor for periodic E-field not implemented") 254 END IF 255 256 fieldpol = dft_control%period_efield%polarisation 257 fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) 258 fieldpol = -fieldpol*dft_control%period_efield%strength 259 hmat = cell%hmat(:, :)/twopi 260 DO idir = 1, 3 261 fpolvec(idir) = fieldpol(1)*hmat(1, idir) + fieldpol(2)*hmat(2, idir) + fieldpol(3)*hmat(3, idir) 262 END DO 263 264 ! nuclear contribution 265 natom = SIZE(particle_set) 266 IF (calculate_forces) THEN 267 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force) 268 ALLOCATE (atom_of_kind(natom)) 269 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 270 END IF 271 zi(:) = CMPLX(1._dp, 0._dp, dp) 272 DO ia = 1, natom 273 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind) 274 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge) 275 ria = particle_set(ia)%r 276 ria = pbc(ria, cell) 277 DO idir = 1, 3 278 kvec(:) = twopi*cell%h_inv(idir, :) 279 dd = SUM(kvec(:)*ria(:)) 280 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge 281 zi(idir) = zi(idir)*zdeta 282 END DO 283 IF (calculate_forces) THEN 284 IF (para_env%mepos == 0) THEN 285 iatom = atom_of_kind(ia) 286 forcea(:) = fieldpol(:)*charge 287 force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) + forcea(:) 288 END IF 289 END IF 290 IF (use_virial) THEN 291 IF (para_env%mepos == 0) & 292 CALL virial_pair_force(virial%pv_virial, 1.0_dp, forcea, ria) 293 END IF 294 END DO 295 qi = AIMAG(LOG(zi)) 296 297 ! check uniform occupation 298 NULLIFY (mos) 299 CALL get_qs_env(qs_env=qs_env, mos=mos) 300 DO ispin = 1, dft_control%nspins 301 CALL get_mo_set(mo_set=mos(ispin)%mo_set, maxocc=occ, uniform_occupation=uniform) 302 IF (.NOT. uniform) THEN 303 CPABORT("Berry phase moments for non uniform MOs' occupation numbers not implemented") 304 END IF 305 END DO 306 307 NULLIFY (mo_derivs) 308 CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs) 309 ! initialize all work matrices needed 310 ALLOCATE (op_fm_set(2, dft_control%nspins)) 311 ALLOCATE (opvec(2, dft_control%nspins)) 312 ALLOCATE (eigrmat(dft_control%nspins)) 313 ALLOCATE (inv_mat(dft_control%nspins)) 314 ALLOCATE (inv_work(2, dft_control%nspins)) 315 ALLOCATE (mo_derivs_tmp(SIZE(mo_derivs))) 316 ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs))) 317 318 ! Allocate temp matrices for the wavefunction derivatives 319 DO ispin = 1, dft_control%nspins 320 NULLIFY (tmp_fm_struct, mo_coeff) 321 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo) 322 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, & 323 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context) 324 CALL cp_fm_create(mo_derivs_tmp(ispin)%matrix, mo_coeff%matrix_struct) 325 CALL cp_fm_create(mo_coeff_tmp(ispin)%matrix, mo_coeff%matrix_struct) 326 CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_derivs_tmp(ispin)%matrix) 327 DO i = 1, SIZE(op_fm_set, 1) 328 CALL cp_fm_create(opvec(i, ispin)%matrix, mo_coeff%matrix_struct) 329 NULLIFY (op_fm_set(i, ispin)%matrix) 330 CALL cp_fm_create(op_fm_set(i, ispin)%matrix, tmp_fm_struct) 331 CALL cp_fm_create(inv_work(i, ispin)%matrix, op_fm_set(i, ispin)%matrix%matrix_struct) 332 END DO 333 CALL cp_cfm_create(eigrmat(ispin)%matrix, op_fm_set(1, ispin)%matrix%matrix_struct) 334 CALL cp_cfm_create(inv_mat(ispin)%matrix, op_fm_set(1, ispin)%matrix%matrix_struct) 335 CALL cp_fm_struct_release(tmp_fm_struct) 336 END DO 337 ! temp matrices for force calculation 338 IF (calculate_forces) THEN 339 NULLIFY (matrix_s) 340 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s) 341 ALLOCATE (tempmat(2, dft_control%nspins)) 342 DO ispin = 1, dft_control%nspins 343 ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix) 344 CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT') 345 CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT') 346 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp) 347 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp) 348 END DO 349 ! integration 350 CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab) 351 ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab)) 352 ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3)) 353 lsab = MAX(ldab, lsab) 354 DO i = 1, 3 355 ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab)) 356 ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab)) 357 END DO 358 END IF 359 360 !Start the MO derivative calculation 361 !loop over all cell vectors 362 DO idir = 1, 3 363 ci(idir) = 0.0_dp 364 zi(idir) = z_zero 365 IF (ABS(fpolvec(idir)) .GT. 1.0E-12_dp) THEN 366 cosmat => efield%cosmat(idir)%matrix 367 sinmat => efield%sinmat(idir)%matrix 368 !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1) 369 !first step S_berry * C and C^T S_berry C 370 DO ispin = 1, dft_control%nspins ! spin 371 IF (mos(ispin)%mo_set%use_mo_coeff_b) THEN 372 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo) 373 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin)%matrix) 374 ELSE 375 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, mo_coeff=mo_coeff_tmp(ispin)%matrix, nmo=nmo) 376 END IF 377 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin)%matrix, opvec(1, ispin)%matrix, ncol=nmo) 378 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin)%matrix, opvec(1, ispin)%matrix, 0.0_dp, & 379 op_fm_set(1, ispin)%matrix) 380 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin)%matrix, opvec(2, ispin)%matrix, ncol=nmo) 381 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin)%matrix, opvec(2, ispin)%matrix, 0.0_dp, & 382 op_fm_set(2, ispin)%matrix) 383 ENDDO 384 !second step invert C^T S_berry C 385 zdet = z_one 386 DO ispin = 1, dft_control%nspins 387 CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin)%matrix, z_one, op_fm_set(1, ispin)%matrix) 388 CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin)%matrix, -gaussi, op_fm_set(2, ispin)%matrix) 389 CALL cp_cfm_set_all(inv_mat(ispin)%matrix, z_zero, z_one) 390 CALL cp_cfm_solve(eigrmat(ispin)%matrix, inv_mat(ispin)%matrix, zdeta) 391 zdet = zdet*zdeta 392 END DO 393 zi(idir) = zdet**occ 394 ci(idir) = AIMAG(LOG(zdet**occ)) 395 396 IF (.NOT. just_energy) THEN 397 !compute the orbital derivative 398 focc = fpolvec(idir) 399 DO ispin = 1, dft_control%nspins 400 inv_work(1, ispin)%matrix%local_data(:, :) = REAL(inv_mat(ispin)%matrix%local_data(:, :), dp) 401 inv_work(2, ispin)%matrix%local_data(:, :) = AIMAG(inv_mat(ispin)%matrix%local_data(:, :)) 402 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) 403 CALL cp_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin)%matrix, inv_work(2, ispin)%matrix, & 404 1.0_dp, mo_derivs_tmp(ispin)%matrix) 405 CALL cp_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin)%matrix, inv_work(1, ispin)%matrix, & 406 1.0_dp, mo_derivs_tmp(ispin)%matrix) 407 END DO 408 END IF 409 410 !compute nuclear forces 411 IF (calculate_forces) THEN 412 nkind = SIZE(qs_kind_set) 413 natom = SIZE(particle_set) 414 kvec(:) = twopi*cell%h_inv(idir, :) 415 416 ! calculate: C [C^T S_berry C]^(-1) C^T 417 ! Store this matrix in DBCSR form (only S overlap blocks) 418 DO ispin = 1, dft_control%nspins 419 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp) 420 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp) 421 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) 422 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin)%matrix, inv_work(1, ispin)%matrix, 0.0_dp, & 423 opvec(1, ispin)%matrix) 424 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin)%matrix, inv_work(2, ispin)%matrix, 0.0_dp, & 425 opvec(2, ispin)%matrix) 426 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, & 427 matrix_v=opvec(1, ispin)%matrix, matrix_g=mo_coeff_tmp(ispin)%matrix, ncol=nmo) 428 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, & 429 matrix_v=opvec(2, ispin)%matrix, matrix_g=mo_coeff_tmp(ispin)%matrix, ncol=nmo) 430 END DO 431 432 ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db) 433 ALLOCATE (basis_set_list(nkind)) 434 DO ikind = 1, nkind 435 qs_kind => qs_kind_set(ikind) 436 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a) 437 IF (ASSOCIATED(basis_set_a)) THEN 438 basis_set_list(ikind)%gto_basis_set => basis_set_a 439 ELSE 440 NULLIFY (basis_set_list(ikind)%gto_basis_set) 441 END IF 442 END DO 443 ! 444 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 445 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 446 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 447 iatom=iatom, jatom=jatom, r=rab) 448 basis_set_a => basis_set_list(ikind)%gto_basis_set 449 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 450 basis_set_b => basis_set_list(jkind)%gto_basis_set 451 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 452 ! basis ikind 453 first_sgfa => basis_set_a%first_sgf 454 la_max => basis_set_a%lmax 455 la_min => basis_set_a%lmin 456 npgfa => basis_set_a%npgf 457 nseta = basis_set_a%nset 458 nsgfa => basis_set_a%nsgf_set 459 rpgfa => basis_set_a%pgf_radius 460 set_radius_a => basis_set_a%set_radius 461 sphi_a => basis_set_a%sphi 462 zeta => basis_set_a%zet 463 ! basis jkind 464 first_sgfb => basis_set_b%first_sgf 465 lb_max => basis_set_b%lmax 466 lb_min => basis_set_b%lmin 467 npgfb => basis_set_b%npgf 468 nsetb = basis_set_b%nset 469 nsgfb => basis_set_b%nsgf_set 470 rpgfb => basis_set_b%pgf_radius 471 set_radius_b => basis_set_b%set_radius 472 sphi_b => basis_set_b%sphi 473 zetb => basis_set_b%zet 474 475 atom_a = atom_of_kind(iatom) 476 atom_b = atom_of_kind(jatom) 477 478 ldsa = SIZE(sphi_a, 1) 479 ldsb = SIZE(sphi_b, 1) 480 ra(:) = pbc(particle_set(iatom)%r(:), cell) 481 rb(:) = ra + rab 482 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 483 484 IF (iatom <= jatom) THEN 485 irow = iatom 486 icol = jatom 487 ELSE 488 irow = jatom 489 icol = iatom 490 END IF 491 492 IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN 493 fab = 1.0_dp*occ 494 ELSE 495 fab = 2.0_dp*occ 496 END IF 497 498 DO i = 1, 3 499 dcost(i, 1)%block = 0.0_dp 500 dsint(i, 1)%block = 0.0_dp 501 dcost(i, 2)%block = 0.0_dp 502 dsint(i, 2)%block = 0.0_dp 503 END DO 504 505 DO iset = 1, nseta 506 ncoa = npgfa(iset)*ncoset(la_max(iset)) 507 sgfa = first_sgfa(1, iset) 508 DO jset = 1, nsetb 509 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 510 ncob = npgfb(jset)*ncoset(lb_max(jset)) 511 sgfb = first_sgfb(1, jset) 512 ! Calculate the primitive integrals (da|b) 513 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 514 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 515 ra, rb, kvec, cosab, sinab, dcosab, dsinab) 516 DO i = 1, 3 517 CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, & 518 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, & 519 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, & 520 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab) 521 END DO 522 ! Calculate the primitive integrals (a|db) 523 CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 524 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 525 rb, ra, kvec, cosab, sinab, dcosab, dsinab) 526 DO i = 1, 3 527 dcosab(1:ncoa, 1:ncob, i) = TRANSPOSE(dcosab(1:ncob, 1:ncoa, i)) 528 dsinab(1:ncoa, 1:ncob, i) = TRANSPOSE(dsinab(1:ncob, 1:ncoa, i)) 529 CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, & 530 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, & 531 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, & 532 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab) 533 END DO 534 END DO 535 END DO 536 forcea = 0.0_dp 537 forceb = 0.0_dp 538 DO ispin = 1, dft_control%nspins 539 NULLIFY (rblock, iblock) 540 CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, & 541 row=irow, col=icol, BLOCK=rblock, found=found) 542 CPASSERT(found) 543 CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, & 544 row=irow, col=icol, BLOCK=iblock, found=found) 545 CPASSERT(found) 546 n1 = SIZE(rblock, 1) 547 n2 = SIZE(rblock, 2) 548 CPASSERT(SIZE(iblock, 1) == n1) 549 CPASSERT(SIZE(iblock, 2) == n2) 550 CPASSERT(lsab >= n1) 551 CPASSERT(lsab >= n2) 552 IF (iatom <= jatom) THEN 553 DO i = 1, 3 554 forcea(i) = forcea(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) & 555 - SUM(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2)) 556 forceb(i) = forceb(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) & 557 - SUM(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2)) 558 END DO 559 ELSE 560 DO i = 1, 3 561 forcea(i) = forcea(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) & 562 - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1)) 563 forceb(i) = forceb(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) & 564 - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1)) 565 END DO 566 END IF 567 END DO 568 force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) - fab*fpolvec(idir)*forcea(1:3) 569 force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) - fab*fpolvec(idir)*forceb(1:3) 570 IF (use_virial) THEN 571 f0 = -fab*fpolvec(idir) 572 CALL virial_pair_force(virial%pv_virial, f0, forcea, ra) 573 CALL virial_pair_force(virial%pv_virial, f0, forceb, rb) 574 END IF 575 576 END DO 577 CALL neighbor_list_iterator_release(nl_iterator) 578 DEALLOCATE (basis_set_list) 579 580 END IF 581 END IF 582 END DO 583 584 ! Energy 585 ener_field = 0.0_dp 586 ti = 0.0_dp 587 DO idir = 1, 3 588 ! make sure the total normalized polarization is within [-1:1] 589 cqi(idir) = qi(idir) + ci(idir) 590 IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi 591 IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi 592 ! now check for log branch 593 IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN 594 ti(idir) = (efield%polarisation(idir) - cqi(idir))/pi 595 DO i = 1, 10 596 cqi(idir) = cqi(idir) + SIGN(1.0_dp, ti(idir))*twopi 597 IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT 598 END DO 599 END IF 600 ener_field = ener_field + fpolvec(idir)*cqi(idir) 601 END DO 602 603 ! update the references 604 IF (calculate_forces) THEN 605 ! check for smoothness of energy surface 606 IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(fpolvec))) THEN 607 CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface") 608 END IF 609 efield%field_energy = ener_field 610 efield%polarisation(:) = cqi(:) 611 END IF 612 energy%efield = ener_field 613 614 IF (.NOT. just_energy) THEN 615 ! Add the result to mo_derivativs 616 DO ispin = 1, dft_control%nspins 617 CALL copy_fm_to_dbcsr(mo_derivs_tmp(ispin)%matrix, mo_derivs(ispin)%matrix) 618 END DO 619 IF (use_virial) THEN 620 ti = 0.0_dp 621 DO i = 1, 3 622 DO j = 1, 3 623 ti(j) = ti(j) + hmat(j, i)*cqi(i) 624 END DO 625 END DO 626 DO i = 1, 3 627 DO j = 1, 3 628 virial%pv_virial(i, j) = virial%pv_virial(i, j) - fieldpol(i)*ti(j) 629 END DO 630 END DO 631 END IF 632 END IF 633 634 DO ispin = 1, dft_control%nspins 635 CALL cp_cfm_release(eigrmat(ispin)%matrix) 636 CALL cp_cfm_release(inv_mat(ispin)%matrix) 637 CALL cp_fm_release(mo_derivs_tmp(ispin)%matrix) 638 CALL cp_fm_release(mo_coeff_tmp(ispin)%matrix) 639 DO i = 1, SIZE(op_fm_set, 1) 640 CALL cp_fm_release(opvec(i, ispin)%matrix) 641 CALL cp_fm_release(op_fm_set(i, ispin)%matrix) 642 CALL cp_fm_release(inv_work(i, ispin)%matrix) 643 END DO 644 END DO 645 DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat) 646 DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp) 647 648 IF (calculate_forces) THEN 649 DO ikind = 1, SIZE(atomic_kind_set) 650 CALL mp_sum(force(ikind)%efield, para_env%group) 651 END DO 652 DEALLOCATE (atom_of_kind) 653 DEALLOCATE (cosab, sinab, work, dcosab, dsinab) 654 DO i = 1, 3 655 DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block) 656 DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block) 657 END DO 658 CALL dbcsr_deallocate_matrix_set(tempmat) 659 END IF 660 CALL timestop(handle) 661 662 END SUBROUTINE qs_efield_derivatives 663 664! ************************************************************************************************** 665!> \brief ... 666!> \param qs_env ... 667!> \param just_energy ... 668!> \param calculate_forces ... 669! ************************************************************************************************** 670 SUBROUTINE qs_dispfield_derivatives(qs_env, just_energy, calculate_forces) 671 TYPE(qs_environment_type), POINTER :: qs_env 672 LOGICAL, INTENT(IN) :: just_energy, calculate_forces 673 674 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_dispfield_derivatives', & 675 routineP = moduleN//':'//routineN 676 677 COMPLEX(dp) :: zdet, zdeta, zi(3) 678 INTEGER :: handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, jatom, jkind, jset, & 679 ldab, ldsa, ldsb, lsab, n1, n2, nao, natom, ncoa, ncob, nkind, nmo, nseta, nsetb, sgfa, & 680 sgfb 681 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 682 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 683 npgfb, nsgfa, nsgfb 684 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 685 LOGICAL :: found, uniform, use_virial 686 REAL(dp) :: charge, ci(3), cqi(3), dab, dd, di(3), & 687 ener_field, fab, fieldpol(3), focc, & 688 hmat(3, 3), occ, omega, qi(3), & 689 rlog(3), zlog(3) 690 REAL(dp), DIMENSION(3) :: dfilter, forcea, forceb, kvec, ra, rab, & 691 rb, ria 692 REAL(dp), DIMENSION(:, :), POINTER :: cosab, iblock, rblock, sinab, work 693 REAL(dp), DIMENSION(:, :, :), POINTER :: dcosab, dsinab, force_tmp 694 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 695 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 696 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 697 TYPE(block_p_type), DIMENSION(3, 2) :: dcost, dsint 698 TYPE(cell_type), POINTER :: cell 699 TYPE(cp_cfm_p_type), DIMENSION(:), POINTER :: eigrmat, inv_mat 700 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_coeff_tmp 701 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: inv_work, mo_derivs_tmp, op_fm_set, opvec 702 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 703 TYPE(cp_fm_type), POINTER :: mo_coeff 704 TYPE(cp_para_env_type), POINTER :: para_env 705 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, mo_derivs 706 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: tempmat 707 TYPE(dbcsr_type), POINTER :: cosmat, mo_coeff_b, sinmat 708 TYPE(dft_control_type), POINTER :: dft_control 709 TYPE(efield_berry_type), POINTER :: efield 710 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 711 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 712 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 713 TYPE(neighbor_list_iterator_p_type), & 714 DIMENSION(:), POINTER :: nl_iterator 715 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 716 POINTER :: sab_orb 717 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 718 TYPE(qs_energy_type), POINTER :: energy 719 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 720 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 721 TYPE(qs_kind_type), POINTER :: qs_kind 722 TYPE(virial_type), POINTER :: virial 723 724 CALL timeset(routineN, handle) 725 726 NULLIFY (dft_control, cell, particle_set) 727 CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, & 728 particle_set=particle_set, virial=virial) 729 NULLIFY (qs_kind_set, efield, para_env, sab_orb) 730 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 731 efield=efield, energy=energy, para_env=para_env, sab_orb=sab_orb) 732 733 ! calculate stress only if forces requested also 734 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 735 use_virial = use_virial .AND. calculate_forces 736 ! disable stress calculation 737 IF (use_virial) THEN 738 CPABORT("Stress tensor for periodic D-field not implemented") 739 END IF 740 741 dfilter(1:3) = dft_control%period_efield%d_filter(1:3) 742 743 fieldpol = dft_control%period_efield%polarisation 744 fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol)) 745 fieldpol = fieldpol*dft_control%period_efield%strength 746 747 omega = cell%deth 748 hmat = cell%hmat(:, :)/(twopi*omega) 749 750 ! nuclear contribution to polarization 751 natom = SIZE(particle_set) 752 IF (calculate_forces) THEN 753 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force) 754 ALLOCATE (atom_of_kind(natom)) 755 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 756 ALLOCATE (force_tmp(natom, 3, 3)) 757 force_tmp = 0.0_dp 758 END IF 759 zi(:) = CMPLX(1._dp, 0._dp, dp) 760 DO ia = 1, natom 761 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind) 762 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge) 763 ria = particle_set(ia)%r 764 ria = pbc(ria, cell) 765 DO idir = 1, 3 766 kvec(:) = twopi*cell%h_inv(idir, :) 767 dd = SUM(kvec(:)*ria(:)) 768 zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)**charge 769 zi(idir) = zi(idir)*zdeta 770 END DO 771 IF (calculate_forces) THEN 772 IF (para_env%mepos == 0) THEN 773 DO i = 1, 3 774 force_tmp(ia, i, i) = force_tmp(ia, i, i) + charge/omega 775 END DO 776 END IF 777 END IF 778 END DO 779 rlog = AIMAG(LOG(zi)) 780 781 ! check uniform occupation 782 NULLIFY (mos) 783 CALL get_qs_env(qs_env=qs_env, mos=mos) 784 DO ispin = 1, dft_control%nspins 785 CALL get_mo_set(mo_set=mos(ispin)%mo_set, maxocc=occ, uniform_occupation=uniform) 786 IF (.NOT. uniform) THEN 787 CPABORT("Berry phase moments for non uniform MO occupation numbers not implemented") 788 END IF 789 END DO 790 791 ! initialize all work matrices needed 792 NULLIFY (mo_derivs) 793 CALL get_qs_env(qs_env=qs_env, mo_derivs=mo_derivs) 794 ALLOCATE (op_fm_set(2, dft_control%nspins)) 795 ALLOCATE (opvec(2, dft_control%nspins)) 796 ALLOCATE (eigrmat(dft_control%nspins)) 797 ALLOCATE (inv_mat(dft_control%nspins)) 798 ALLOCATE (inv_work(2, dft_control%nspins)) 799 ALLOCATE (mo_derivs_tmp(3, SIZE(mo_derivs))) 800 ALLOCATE (mo_coeff_tmp(SIZE(mo_derivs))) 801 802 ! Allocate temp matrices for the wavefunction derivatives 803 DO ispin = 1, dft_control%nspins 804 NULLIFY (tmp_fm_struct, mo_coeff) 805 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo) 806 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, & 807 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context) 808 CALL cp_fm_create(mo_coeff_tmp(ispin)%matrix, mo_coeff%matrix_struct) 809 DO i = 1, 3 810 CALL cp_fm_create(mo_derivs_tmp(i, ispin)%matrix, mo_coeff%matrix_struct) 811 CALL cp_fm_set_all(matrix=mo_derivs_tmp(i, ispin)%matrix, alpha=0.0_dp) 812 END DO 813 DO i = 1, SIZE(op_fm_set, 1) 814 CALL cp_fm_create(opvec(i, ispin)%matrix, mo_coeff%matrix_struct) 815 NULLIFY (op_fm_set(i, ispin)%matrix) 816 CALL cp_fm_create(op_fm_set(i, ispin)%matrix, tmp_fm_struct) 817 CALL cp_fm_create(inv_work(i, ispin)%matrix, op_fm_set(i, ispin)%matrix%matrix_struct) 818 END DO 819 CALL cp_cfm_create(eigrmat(ispin)%matrix, op_fm_set(1, ispin)%matrix%matrix_struct) 820 CALL cp_cfm_create(inv_mat(ispin)%matrix, op_fm_set(1, ispin)%matrix%matrix_struct) 821 CALL cp_fm_struct_release(tmp_fm_struct) 822 END DO 823 ! temp matrices for force calculation 824 IF (calculate_forces) THEN 825 NULLIFY (matrix_s) 826 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s) 827 ALLOCATE (tempmat(2, dft_control%nspins)) 828 DO ispin = 1, dft_control%nspins 829 ALLOCATE (tempmat(1, ispin)%matrix, tempmat(2, ispin)%matrix) 830 CALL dbcsr_copy(tempmat(1, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT') 831 CALL dbcsr_copy(tempmat(2, ispin)%matrix, matrix_s(1)%matrix, 'TEMPMAT') 832 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp) 833 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp) 834 END DO 835 ! integration 836 CALL get_qs_kind_set(qs_kind_set, maxco=ldab, maxsgf=lsab) 837 ALLOCATE (cosab(ldab, ldab), sinab(ldab, ldab), work(ldab, ldab)) 838 ALLOCATE (dcosab(ldab, ldab, 3), dsinab(ldab, ldab, 3)) 839 lsab = MAX(lsab, ldab) 840 DO i = 1, 3 841 ALLOCATE (dcost(i, 1)%block(lsab, lsab), dsint(i, 1)%block(lsab, lsab)) 842 ALLOCATE (dcost(i, 2)%block(lsab, lsab), dsint(i, 2)%block(lsab, lsab)) 843 END DO 844 END IF 845 846 !Start the MO derivative calculation 847 !loop over all cell vectors 848 DO idir = 1, 3 849 zi(idir) = z_zero 850 cosmat => efield%cosmat(idir)%matrix 851 sinmat => efield%sinmat(idir)%matrix 852 !evaluate the expression needed for the derivative (S_berry * C and [C^T S_berry C]^-1) 853 !first step S_berry * C and C^T S_berry C 854 DO ispin = 1, dft_control%nspins ! spin 855 IF (mos(ispin)%mo_set%use_mo_coeff_b) THEN 856 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, mo_coeff_b=mo_coeff_b, nmo=nmo) 857 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff_tmp(ispin)%matrix) 858 ELSE 859 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, mo_coeff=mo_coeff_tmp(ispin)%matrix, nmo=nmo) 860 END IF 861 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff_tmp(ispin)%matrix, opvec(1, ispin)%matrix, ncol=nmo) 862 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin)%matrix, opvec(1, ispin)%matrix, 0.0_dp, & 863 op_fm_set(1, ispin)%matrix) 864 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff_tmp(ispin)%matrix, opvec(2, ispin)%matrix, ncol=nmo) 865 CALL cp_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff_tmp(ispin)%matrix, opvec(2, ispin)%matrix, 0.0_dp, & 866 op_fm_set(2, ispin)%matrix) 867 ENDDO 868 !second step invert C^T S_berry C 869 zdet = z_one 870 DO ispin = 1, dft_control%nspins 871 CALL cp_cfm_scale_and_add_fm(z_zero, eigrmat(ispin)%matrix, z_one, op_fm_set(1, ispin)%matrix) 872 CALL cp_cfm_scale_and_add_fm(z_one, eigrmat(ispin)%matrix, -gaussi, op_fm_set(2, ispin)%matrix) 873 CALL cp_cfm_set_all(inv_mat(ispin)%matrix, z_zero, z_one) 874 CALL cp_cfm_solve(eigrmat(ispin)%matrix, inv_mat(ispin)%matrix, zdeta) 875 zdet = zdet*zdeta 876 END DO 877 zi(idir) = zdet**occ 878 zlog(idir) = AIMAG(LOG(zi(idir))) 879 880 IF (.NOT. just_energy) THEN 881 !compute the orbital derivative 882 DO ispin = 1, dft_control%nspins 883 inv_work(1, ispin)%matrix%local_data(:, :) = REAL(inv_mat(ispin)%matrix%local_data(:, :), dp) 884 inv_work(2, ispin)%matrix%local_data(:, :) = AIMAG(inv_mat(ispin)%matrix%local_data(:, :)) 885 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) 886 DO i = 1, 3 887 focc = hmat(idir, i) 888 CALL cp_gemm("N", "N", nao, nmo, nmo, focc, opvec(1, ispin)%matrix, inv_work(2, ispin)%matrix, & 889 1.0_dp, mo_derivs_tmp(idir, ispin)%matrix) 890 CALL cp_gemm("N", "N", nao, nmo, nmo, -focc, opvec(2, ispin)%matrix, inv_work(1, ispin)%matrix, & 891 1.0_dp, mo_derivs_tmp(idir, ispin)%matrix) 892 END DO 893 END DO 894 END IF 895 896 !compute nuclear forces 897 IF (calculate_forces) THEN 898 nkind = SIZE(qs_kind_set) 899 natom = SIZE(particle_set) 900 kvec(:) = twopi*cell%h_inv(idir, :) 901 902 ! calculate: C [C^T S_berry C]^(-1) C^T 903 ! Store this matrix in DBCSR form (only S overlap blocks) 904 DO ispin = 1, dft_control%nspins 905 CALL dbcsr_set(tempmat(1, ispin)%matrix, 0.0_dp) 906 CALL dbcsr_set(tempmat(2, ispin)%matrix, 0.0_dp) 907 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) 908 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin)%matrix, inv_work(1, ispin)%matrix, 0.0_dp, & 909 opvec(1, ispin)%matrix) 910 CALL cp_gemm("N", "N", nao, nmo, nmo, 1.0_dp, mo_coeff_tmp(ispin)%matrix, inv_work(2, ispin)%matrix, 0.0_dp, & 911 opvec(2, ispin)%matrix) 912 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(1, ispin)%matrix, & 913 matrix_v=opvec(1, ispin)%matrix, matrix_g=mo_coeff_tmp(ispin)%matrix, ncol=nmo) 914 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=tempmat(2, ispin)%matrix, & 915 matrix_v=opvec(2, ispin)%matrix, matrix_g=mo_coeff_tmp(ispin)%matrix, ncol=nmo) 916 END DO 917 918 ! Calculation of derivative integrals (da|eikr|b) and (a|eikr|db) 919 ALLOCATE (basis_set_list(nkind)) 920 DO ikind = 1, nkind 921 qs_kind => qs_kind_set(ikind) 922 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a) 923 IF (ASSOCIATED(basis_set_a)) THEN 924 basis_set_list(ikind)%gto_basis_set => basis_set_a 925 ELSE 926 NULLIFY (basis_set_list(ikind)%gto_basis_set) 927 END IF 928 END DO 929 ! 930 CALL neighbor_list_iterator_create(nl_iterator, sab_orb) 931 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 932 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, & 933 iatom=iatom, jatom=jatom, r=rab) 934 basis_set_a => basis_set_list(ikind)%gto_basis_set 935 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 936 basis_set_b => basis_set_list(jkind)%gto_basis_set 937 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 938 ! basis ikind 939 first_sgfa => basis_set_a%first_sgf 940 la_max => basis_set_a%lmax 941 la_min => basis_set_a%lmin 942 npgfa => basis_set_a%npgf 943 nseta = basis_set_a%nset 944 nsgfa => basis_set_a%nsgf_set 945 rpgfa => basis_set_a%pgf_radius 946 set_radius_a => basis_set_a%set_radius 947 sphi_a => basis_set_a%sphi 948 zeta => basis_set_a%zet 949 ! basis jkind 950 first_sgfb => basis_set_b%first_sgf 951 lb_max => basis_set_b%lmax 952 lb_min => basis_set_b%lmin 953 npgfb => basis_set_b%npgf 954 nsetb = basis_set_b%nset 955 nsgfb => basis_set_b%nsgf_set 956 rpgfb => basis_set_b%pgf_radius 957 set_radius_b => basis_set_b%set_radius 958 sphi_b => basis_set_b%sphi 959 zetb => basis_set_b%zet 960 961 ldsa = SIZE(sphi_a, 1) 962 ldsb = SIZE(sphi_b, 1) 963 ra(:) = pbc(particle_set(iatom)%r(:), cell) 964 rb(:) = ra + rab 965 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 966 967 IF (iatom <= jatom) THEN 968 irow = iatom 969 icol = jatom 970 ELSE 971 irow = jatom 972 icol = iatom 973 END IF 974 975 IF (iatom == jatom .AND. dab < 1.e-10_dp) THEN 976 fab = 1.0_dp*occ 977 ELSE 978 fab = 2.0_dp*occ 979 END IF 980 981 DO i = 1, 3 982 dcost(i, 1)%block = 0.0_dp 983 dsint(i, 1)%block = 0.0_dp 984 dcost(i, 2)%block = 0.0_dp 985 dsint(i, 2)%block = 0.0_dp 986 END DO 987 988 DO iset = 1, nseta 989 ncoa = npgfa(iset)*ncoset(la_max(iset)) 990 sgfa = first_sgfa(1, iset) 991 DO jset = 1, nsetb 992 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 993 ncob = npgfb(jset)*ncoset(lb_max(jset)) 994 sgfb = first_sgfb(1, jset) 995 ! Calculate the primitive integrals (da|b) 996 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 997 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 998 ra, rb, kvec, cosab, sinab, dcosab, dsinab) 999 DO i = 1, 3 1000 CALL contract_all(dcost(i, 1)%block, dsint(i, 1)%block, & 1001 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, & 1002 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, & 1003 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab) 1004 END DO 1005 ! Calculate the primitive integrals (a|db) 1006 CALL cossin(lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), & 1007 la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), & 1008 rb, ra, kvec, cosab, sinab, dcosab, dsinab) 1009 DO i = 1, 3 1010 dcosab(1:ncoa, 1:ncob, i) = TRANSPOSE(dcosab(1:ncob, 1:ncoa, i)) 1011 dsinab(1:ncoa, 1:ncob, i) = TRANSPOSE(dsinab(1:ncob, 1:ncoa, i)) 1012 CALL contract_all(dcost(i, 2)%block, dsint(i, 2)%block, & 1013 ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, & 1014 ncob, nsgfb(jset), sgfb, sphi_b, ldsb, & 1015 dcosab(:, :, i), dsinab(:, :, i), ldab, work, ldab) 1016 END DO 1017 END DO 1018 END DO 1019 forcea = 0.0_dp 1020 forceb = 0.0_dp 1021 DO ispin = 1, dft_control%nspins 1022 NULLIFY (rblock, iblock) 1023 CALL dbcsr_get_block_p(matrix=tempmat(1, ispin)%matrix, & 1024 row=irow, col=icol, BLOCK=rblock, found=found) 1025 CPASSERT(found) 1026 CALL dbcsr_get_block_p(matrix=tempmat(2, ispin)%matrix, & 1027 row=irow, col=icol, BLOCK=iblock, found=found) 1028 CPASSERT(found) 1029 n1 = SIZE(rblock, 1) 1030 n2 = SIZE(rblock, 2) 1031 CPASSERT(SIZE(iblock, 1) == n1) 1032 CPASSERT(SIZE(iblock, 2) == n2) 1033 CPASSERT(lsab >= n1) 1034 CPASSERT(lsab >= n2) 1035 IF (iatom <= jatom) THEN 1036 DO i = 1, 3 1037 forcea(i) = forcea(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 1)%block(1:n1, 1:n2)) & 1038 - SUM(iblock(1:n1, 1:n2)*dcost(i, 1)%block(1:n1, 1:n2)) 1039 forceb(i) = forceb(i) + SUM(rblock(1:n1, 1:n2)*dsint(i, 2)%block(1:n1, 1:n2)) & 1040 - SUM(iblock(1:n1, 1:n2)*dcost(i, 2)%block(1:n1, 1:n2)) 1041 END DO 1042 ELSE 1043 DO i = 1, 3 1044 forcea(i) = forcea(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 1)%block(1:n2, 1:n1)) & 1045 - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 1)%block(1:n2, 1:n1)) 1046 forceb(i) = forceb(i) + SUM(TRANSPOSE(rblock(1:n1, 1:n2))*dsint(i, 2)%block(1:n2, 1:n1)) & 1047 - SUM(TRANSPOSE(iblock(1:n1, 1:n2))*dcost(i, 2)%block(1:n2, 1:n1)) 1048 END DO 1049 END IF 1050 END DO 1051 DO i = 1, 3 1052 force_tmp(iatom, :, i) = force_tmp(iatom, :, i) - fab*hmat(i, idir)*forcea(:) 1053 force_tmp(jatom, :, i) = force_tmp(jatom, :, i) - fab*hmat(i, idir)*forceb(:) 1054 END DO 1055 END DO 1056 CALL neighbor_list_iterator_release(nl_iterator) 1057 DEALLOCATE (basis_set_list) 1058 END IF 1059 END DO 1060 1061 ! make sure the total normalized polarization is within [-1:1] 1062 DO idir = 1, 3 1063 cqi(idir) = rlog(idir) + zlog(idir) 1064 IF (cqi(idir) > pi) cqi(idir) = cqi(idir) - twopi 1065 IF (cqi(idir) < -pi) cqi(idir) = cqi(idir) + twopi 1066 ! now check for log branch 1067 IF (calculate_forces) THEN 1068 IF (ABS(efield%polarisation(idir) - cqi(idir)) > pi) THEN 1069 di(idir) = (efield%polarisation(idir) - cqi(idir))/pi 1070 DO i = 1, 10 1071 cqi(idir) = cqi(idir) + SIGN(1.0_dp, di(idir))*twopi 1072 IF (ABS(efield%polarisation(idir) - cqi(idir)) < pi) EXIT 1073 END DO 1074 END IF 1075 END IF 1076 END DO 1077 DO idir = 1, 3 1078 qi(idir) = 0.0_dp 1079 ci(idir) = 0.0_dp 1080 DO i = 1, 3 1081 ci(idir) = ci(idir) + hmat(idir, i)*cqi(i) 1082 END DO 1083 END DO 1084 1085 ! update the references 1086 IF (calculate_forces) THEN 1087 ener_field = SUM(ci) 1088 ! check for smoothness of energy surface 1089 IF (ABS(efield%field_energy - ener_field) > pi*ABS(SUM(hmat))) THEN 1090 CPWARN("Large change of e-field energy detected. Correct for non-smooth energy surface") 1091 END IF 1092 efield%field_energy = ener_field 1093 efield%polarisation(:) = cqi(:) 1094 END IF 1095 1096 ! Energy 1097 ener_field = 0.0_dp 1098 DO i = 1, 3 1099 ener_field = ener_field + dfilter(i)*(fieldpol(i) - 2._dp*twopi*ci(i))**2 1100 END DO 1101 energy%efield = 0.25_dp*omega/twopi*ener_field 1102 1103 IF (.NOT. just_energy) THEN 1104 DO i = 1, 3 1105 di(i) = -omega*(fieldpol(i) - 2._dp*twopi*ci(i))*dfilter(i) 1106 END DO 1107 ! Add the result to mo_derivativs 1108 DO ispin = 1, dft_control%nspins 1109 CALL copy_dbcsr_to_fm(mo_derivs(ispin)%matrix, mo_coeff_tmp(ispin)%matrix) 1110 DO idir = 1, 3 1111 CALL cp_fm_scale_and_add(1.0_dp, mo_coeff_tmp(ispin)%matrix, di(idir), & 1112 mo_derivs_tmp(idir, ispin)%matrix) 1113 END DO 1114 END DO 1115 DO ispin = 1, dft_control%nspins 1116 CALL copy_fm_to_dbcsr(mo_coeff_tmp(ispin)%matrix, mo_derivs(ispin)%matrix) 1117 END DO 1118 END IF 1119 1120 IF (calculate_forces) THEN 1121 DO i = 1, 3 1122 DO ia = 1, natom 1123 CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind) 1124 iatom = atom_of_kind(ia) 1125 force(ikind)%efield(1:3, iatom) = force(ikind)%efield(1:3, iatom) + di(i)*force_tmp(ia, 1:3, i) 1126 END DO 1127 END DO 1128 END IF 1129 1130 DO ispin = 1, dft_control%nspins 1131 CALL cp_cfm_release(eigrmat(ispin)%matrix) 1132 CALL cp_cfm_release(inv_mat(ispin)%matrix) 1133 CALL cp_fm_release(mo_coeff_tmp(ispin)%matrix) 1134 DO i = 1, 3 1135 CALL cp_fm_release(mo_derivs_tmp(i, ispin)%matrix) 1136 END DO 1137 DO i = 1, SIZE(op_fm_set, 1) 1138 CALL cp_fm_release(opvec(i, ispin)%matrix) 1139 CALL cp_fm_release(op_fm_set(i, ispin)%matrix) 1140 CALL cp_fm_release(inv_work(i, ispin)%matrix) 1141 END DO 1142 END DO 1143 DEALLOCATE (inv_mat, inv_work, op_fm_set, opvec, eigrmat) 1144 DEALLOCATE (mo_coeff_tmp, mo_derivs_tmp) 1145 1146 IF (calculate_forces) THEN 1147 DO ikind = 1, SIZE(atomic_kind_set) 1148 CALL mp_sum(force(ikind)%efield, para_env%group) 1149 END DO 1150 DEALLOCATE (atom_of_kind) 1151 DEALLOCATE (force_tmp) 1152 DEALLOCATE (cosab, sinab, work, dcosab, dsinab) 1153 DO i = 1, 3 1154 DEALLOCATE (dcost(i, 1)%block, dsint(i, 1)%block) 1155 DEALLOCATE (dcost(i, 2)%block, dsint(i, 2)%block) 1156 END DO 1157 CALL dbcsr_deallocate_matrix_set(tempmat) 1158 END IF 1159 CALL timestop(handle) 1160 1161 END SUBROUTINE qs_dispfield_derivatives 1162 1163! ************************************************************************************************** 1164!> \brief ... 1165!> \param cos_block ... 1166!> \param sin_block ... 1167!> \param ncoa ... 1168!> \param nsgfa ... 1169!> \param sgfa ... 1170!> \param sphi_a ... 1171!> \param ldsa ... 1172!> \param ncob ... 1173!> \param nsgfb ... 1174!> \param sgfb ... 1175!> \param sphi_b ... 1176!> \param ldsb ... 1177!> \param cosab ... 1178!> \param sinab ... 1179!> \param ldab ... 1180!> \param work ... 1181!> \param ldwork ... 1182! ************************************************************************************************** 1183 SUBROUTINE contract_all(cos_block, sin_block, & 1184 ncoa, nsgfa, sgfa, sphi_a, ldsa, & 1185 ncob, nsgfb, sgfb, sphi_b, ldsb, & 1186 cosab, sinab, ldab, work, ldwork) 1187 1188 REAL(dp), DIMENSION(:, :), POINTER :: cos_block, sin_block 1189 INTEGER, INTENT(IN) :: ncoa, nsgfa, sgfa 1190 REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_a 1191 INTEGER, INTENT(IN) :: ldsa, ncob, nsgfb, sgfb 1192 REAL(dp), DIMENSION(:, :), INTENT(IN) :: sphi_b 1193 INTEGER, INTENT(IN) :: ldsb 1194 REAL(dp), DIMENSION(:, :), INTENT(IN) :: cosab, sinab 1195 INTEGER, INTENT(IN) :: ldab 1196 REAL(dp), DIMENSION(:, :) :: work 1197 INTEGER, INTENT(IN) :: ldwork 1198 1199! Calculate cosine 1200 1201 CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, cosab(1, 1), ldab, & 1202 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork) 1203 1204 CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, & 1205 work(1, 1), ldwork, 1.0_dp, cos_block(sgfa, sgfb), SIZE(cos_block, 1)) 1206 1207 ! Calculate sine 1208 CALL dgemm("N", "N", ncoa, nsgfb, ncob, 1.0_dp, sinab(1, 1), ldab, & 1209 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork) 1210 1211 CALL dgemm("T", "N", nsgfa, nsgfb, ncoa, 1.0_dp, sphi_a(1, sgfa), ldsa, & 1212 work(1, 1), ldwork, 1.0_dp, sin_block(sgfa, sgfb), SIZE(sin_block, 1)) 1213 1214 END SUBROUTINE contract_all 1215 1216END MODULE qs_efield_berry 1217