1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of the energies concerning the core charge distribution 8!> \par History 9!> - Full refactoring of calculate_ecore and calculate_ecore_overlap (jhu) 10!> \author Matthias Krack (27.04.2001) 11! ************************************************************************************************** 12MODULE qs_core_energies 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind,& 15 get_atomic_kind_set 16 USE atprop_types, ONLY: atprop_array_init,& 17 atprop_type 18 USE cp_para_types, ONLY: cp_para_env_type 19 USE dbcsr_api, ONLY: dbcsr_dot,& 20 dbcsr_p_type,& 21 dbcsr_type 22 USE distribution_1d_types, ONLY: distribution_1d_type 23 USE kinds, ONLY: dp 24 USE mathconstants, ONLY: oorootpi,& 25 twopi 26 USE message_passing, ONLY: mp_sum 27 USE particle_types, ONLY: particle_type 28 USE qs_energy_types, ONLY: qs_energy_type 29 USE qs_environment_types, ONLY: get_qs_env,& 30 qs_environment_type 31 USE qs_force_types, ONLY: qs_force_type 32 USE qs_kind_types, ONLY: get_qs_kind,& 33 qs_kind_type 34 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 35 neighbor_list_iterate,& 36 neighbor_list_iterator_create,& 37 neighbor_list_iterator_p_type,& 38 neighbor_list_iterator_release,& 39 neighbor_list_set_p_type 40 USE virial_methods, ONLY: virial_pair_force 41 USE virial_types, ONLY: virial_type 42#include "./base/base_uses.f90" 43 44 IMPLICIT NONE 45 46 PRIVATE 47 48! *** Global parameters *** 49 50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_energies' 51 52 PUBLIC :: calculate_ptrace, & 53 calculate_ecore_overlap, & 54 calculate_ecore_self 55 56 INTERFACE calculate_ptrace 57 MODULE PROCEDURE calculate_ptrace_1, calculate_ptrace_gamma, calculate_ptrace_kp 58 END INTERFACE 59 60! ************************************************************************************************** 61 62CONTAINS 63 64! ************************************************************************************************** 65!> \brief Calculate the trace of a operator matrix with the density matrix. 66!> Sum over all spin components (in P, no spin in H) 67!> \param hmat ... 68!> \param pmat ... 69!> \param ecore ... 70!> \param nspin ... 71!> \date 29.07.2014 72!> \par History 73!> - none 74!> \author JGH 75!> \version 1.0 76! ************************************************************************************************** 77 SUBROUTINE calculate_ptrace_gamma(hmat, pmat, ecore, nspin) 78 79 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: hmat, pmat 80 REAL(KIND=dp), INTENT(OUT) :: ecore 81 INTEGER, INTENT(IN) :: nspin 82 83 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ptrace_gamma', & 84 routineP = moduleN//':'//routineN 85 86 INTEGER :: handle, ispin 87 REAL(KIND=dp) :: etr 88 89 CALL timeset(routineN, handle) 90 91 ecore = 0.0_dp 92 DO ispin = 1, nspin 93 etr = 0.0_dp 94 CALL dbcsr_dot(hmat(1)%matrix, pmat(ispin)%matrix, etr) 95 ecore = ecore + etr 96 END DO 97 98 CALL timestop(handle) 99 100 END SUBROUTINE calculate_ptrace_gamma 101 102! ************************************************************************************************** 103!> \brief Calculate the trace of a operator matrix with the density matrix. 104!> Sum over all spin components (in P, no spin in H) and the real space 105!> coordinates 106!> \param hmat H matrix 107!> \param pmat P matrices 108!> \param ecore Tr(HP) output 109!> \param nspin Number of P matrices 110!> \date 29.07.2014 111!> \author JGH 112!> \version 1.0 113! ************************************************************************************************** 114 SUBROUTINE calculate_ptrace_kp(hmat, pmat, ecore, nspin) 115 116 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: hmat, pmat 117 REAL(KIND=dp), INTENT(OUT) :: ecore 118 INTEGER, INTENT(IN) :: nspin 119 120 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ptrace_kp', & 121 routineP = moduleN//':'//routineN 122 123 INTEGER :: handle, ic, ispin, nc 124 REAL(KIND=dp) :: etr 125 126 CALL timeset(routineN, handle) 127 128 nc = SIZE(pmat, 2) 129 130 ecore = 0.0_dp 131 DO ispin = 1, nspin 132 DO ic = 1, nc 133 etr = 0.0_dp 134 CALL dbcsr_dot(hmat(1, ic)%matrix, pmat(ispin, ic)%matrix, etr) 135 ecore = ecore + etr 136 END DO 137 END DO 138 139 CALL timestop(handle) 140 141 END SUBROUTINE calculate_ptrace_kp 142 143! ************************************************************************************************** 144!> \brief Calculate the core Hamiltonian energy which includes the kinetic 145!> and the potential energy of the electrons. It is assumed, that 146!> the core Hamiltonian matrix h and the density matrix p have the 147!> same sparse matrix structure (same atomic blocks and block 148!> ordering) 149!> \param h ... 150!> \param p ... 151!> \param ecore ... 152!> \date 03.05.2001 153!> \par History 154!> - simplified taking advantage of new non-redundant matrix 155!> structure (27.06.2003,MK) 156!> - simplified using DBCSR trace function (21.07.2010, jhu) 157!> \author MK 158!> \version 1.0 159! ************************************************************************************************** 160 SUBROUTINE calculate_ptrace_1(h, p, ecore) 161 162 TYPE(dbcsr_type), POINTER :: h, p 163 REAL(KIND=dp), INTENT(OUT) :: ecore 164 165 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ptrace_1', & 166 routineP = moduleN//':'//routineN 167 168 INTEGER :: handle 169 170 CALL timeset(routineN, handle) 171 172 ecore = 0.0_dp 173 CALL dbcsr_dot(h, p, ecore) 174 175 CALL timestop(handle) 176 177 END SUBROUTINE calculate_ptrace_1 178 179! ************************************************************************************************** 180!> \brief Calculate the overlap energy of the core charge distribution. 181!> \param qs_env ... 182!> \param para_env ... 183!> \param calculate_forces ... 184!> \param molecular ... 185!> \param E_overlap_core ... 186!> \date 30.04.2001 187!> \par History 188!> - Force calculation added (03.06.2002,MK) 189!> - Parallelized using a list of local atoms for rows and 190!> columns (19.07.2003,MK) 191!> - Use precomputed neighborlists (sab_core) and nl iterator (28.07.2010,jhu) 192!> \author MK 193!> \version 1.0 194! ************************************************************************************************** 195 SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, & 196 E_overlap_core) 197 TYPE(qs_environment_type), POINTER :: qs_env 198 TYPE(cp_para_env_type), POINTER :: para_env 199 LOGICAL, INTENT(IN) :: calculate_forces 200 LOGICAL, INTENT(IN), OPTIONAL :: molecular 201 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: E_overlap_core 202 203 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_overlap', & 204 routineP = moduleN//':'//routineN 205 206 INTEGER :: atom_a, atom_b, group, handle, iatom, & 207 ikind, jatom, jkind, natom, nkind 208 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 209 LOGICAL :: atenergy, only_molecule, use_virial 210 REAL(KIND=dp) :: aab, dab, eab, ecore_overlap, f, fab, & 211 rab2, rootaab, zab 212 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, radius, zeff 213 REAL(KIND=dp), DIMENSION(3) :: deab, rab 214 REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc 215 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 216 TYPE(atprop_type), POINTER :: atprop 217 TYPE(neighbor_list_iterator_p_type), & 218 DIMENSION(:), POINTER :: nl_iterator 219 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 220 POINTER :: sab_core 221 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 222 TYPE(qs_energy_type), POINTER :: energy 223 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 224 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 225 TYPE(virial_type), POINTER :: virial 226 227 CALL timeset(routineN, handle) 228 229 NULLIFY (atomic_kind_set) 230 NULLIFY (qs_kind_set) 231 NULLIFY (energy) 232 NULLIFY (atprop) 233 NULLIFY (force) 234 NULLIFY (particle_set) 235 236 group = para_env%group 237 238 only_molecule = .FALSE. 239 IF (PRESENT(molecular)) only_molecule = molecular 240 241 CALL get_qs_env(qs_env=qs_env, & 242 atomic_kind_set=atomic_kind_set, & 243 qs_kind_set=qs_kind_set, & 244 particle_set=particle_set, & 245 energy=energy, & 246 force=force, & 247 sab_core=sab_core, & 248 atprop=atprop, & 249 virial=virial) 250 251 ! Allocate work storage 252 nkind = SIZE(atomic_kind_set) 253 natom = SIZE(particle_set) 254 255 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 256 257 ALLOCATE (alpha(nkind), radius(nkind), zeff(nkind)) 258 alpha(:) = 0.0_dp 259 radius(:) = 0.0_dp 260 zeff(:) = 0.0_dp 261 262 IF (calculate_forces) THEN 263 ALLOCATE (atom_of_kind(natom)) 264 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 265 END IF 266 267 atenergy = .FALSE. 268 IF (ASSOCIATED(atprop)) THEN 269 IF (atprop%energy) THEN 270 atenergy = .TRUE. 271 CALL atprop_array_init(atprop%atecc, natom) 272 END IF 273 END IF 274 275 DO ikind = 1, nkind 276 CALL get_qs_kind(qs_kind_set(ikind), & 277 alpha_core_charge=alpha(ikind), & 278 core_charge_radius=radius(ikind), & 279 zeff=zeff(ikind)) 280 END DO 281 282 ecore_overlap = 0.0_dp 283 pv_loc = 0.0_dp 284 285 CALL neighbor_list_iterator_create(nl_iterator, sab_core) 286 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 287 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab) 288 zab = zeff(ikind)*zeff(jkind) 289 aab = alpha(ikind)*alpha(jkind)/(alpha(ikind) + alpha(jkind)) 290 rootaab = SQRT(aab) 291 fab = 2.0_dp*oorootpi*zab*rootaab 292 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 293 IF (rab2 > 1.e-8_dp) THEN 294 IF (ikind == jkind .AND. iatom == jatom) THEN 295 f = 0.5_dp 296 ELSE 297 f = 1.0_dp 298 END IF 299 dab = SQRT(rab2) 300 eab = zab*erfc(rootaab*dab)/dab 301 ecore_overlap = ecore_overlap + f*eab 302 IF (atenergy) THEN 303 atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*f*eab 304 atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*f*eab 305 END IF 306 IF (calculate_forces) THEN 307 deab(:) = rab(:)*f*(eab + fab*EXP(-aab*rab2))/rab2 308 atom_a = atom_of_kind(iatom) 309 atom_b = atom_of_kind(jatom) 310 force(ikind)%core_overlap(:, atom_a) = force(ikind)%core_overlap(:, atom_a) + deab(:) 311 force(jkind)%core_overlap(:, atom_b) = force(jkind)%core_overlap(:, atom_b) - deab(:) 312 IF (use_virial) THEN 313 CALL virial_pair_force(pv_loc, 1._dp, deab, rab) 314 END IF 315 END IF 316 END IF 317 END DO 318 CALL neighbor_list_iterator_release(nl_iterator) 319 320 DEALLOCATE (alpha, radius, zeff) 321 IF (calculate_forces) THEN 322 DEALLOCATE (atom_of_kind) 323 END IF 324 IF (calculate_forces .AND. use_virial) THEN 325 virial%pv_virial = virial%pv_virial + pv_loc 326 virial%pv_hartree = virial%pv_hartree + pv_loc 327 END IF 328 329 CALL mp_sum(ecore_overlap, group) 330 331 energy%core_overlap = ecore_overlap 332 333 IF (PRESENT(E_overlap_core)) THEN 334 E_overlap_core = energy%core_overlap 335 END IF 336 337 CALL timestop(handle) 338 339 END SUBROUTINE calculate_ecore_overlap 340 341! ************************************************************************************************** 342!> \brief Calculate the self energy of the core charge distribution. 343!> \param qs_env ... 344!> \param E_self_core ... 345!> \date 27.04.2001 346!> \author MK 347!> \version 1.0 348! ************************************************************************************************** 349 SUBROUTINE calculate_ecore_self(qs_env, E_self_core) 350 TYPE(qs_environment_type), POINTER :: qs_env 351 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: E_self_core 352 353 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_self', & 354 routineP = moduleN//':'//routineN 355 356 INTEGER :: handle, iatom, ikind, iparticle_local, & 357 natom, nparticle_local 358 REAL(KIND=dp) :: alpha_core_charge, ecore_self, es, zeff 359 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 360 TYPE(atprop_type), POINTER :: atprop 361 TYPE(distribution_1d_type), POINTER :: local_particles 362 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 363 TYPE(qs_energy_type), POINTER :: energy 364 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 365 366! ------------------------------------------------------------------------- 367 368 NULLIFY (atprop) 369 CALL timeset(routineN, handle) 370 371 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, & 372 qs_kind_set=qs_kind_set, energy=energy, atprop=atprop) 373 374 ecore_self = 0.0_dp 375 376 DO ikind = 1, SIZE(atomic_kind_set) 377 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom) 378 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge) 379 ecore_self = ecore_self - REAL(natom, dp)*zeff**2*SQRT(alpha_core_charge) 380 END DO 381 382 energy%core_self = ecore_self/SQRT(twopi) 383 IF (PRESENT(E_self_core)) THEN 384 E_self_core = energy%core_self 385 END IF 386 387 IF (ASSOCIATED(atprop)) THEN 388 IF (atprop%energy) THEN 389 ! atomic energy 390 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, & 391 local_particles=local_particles) 392 natom = SIZE(particle_set) 393 CALL atprop_array_init(atprop%ateself, natom) 394 395 DO ikind = 1, SIZE(atomic_kind_set) 396 nparticle_local = local_particles%n_el(ikind) 397 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge) 398 es = zeff**2*SQRT(alpha_core_charge)/SQRT(twopi) 399 DO iparticle_local = 1, nparticle_local 400 iatom = local_particles%list(ikind)%array(iparticle_local) 401 atprop%ateself(iatom) = atprop%ateself(iatom) - es 402 END DO 403 END DO 404 END IF 405 END IF 406 407 CALL timestop(handle) 408 409 END SUBROUTINE calculate_ecore_self 410 411END MODULE qs_core_energies 412