Line data    Source code

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