LCOV - code coverage report
Current view: top level - src - qs_core_energies.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 100.0 % 174 174
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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 cp_dbcsr_api,                    ONLY: dbcsr_p_type,&
      23              :                                               dbcsr_type
      24              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot
      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          132 :    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          132 :       CALL timeset(routineN, handle)
      93              : 
      94          132 :       ecore = 0.0_dp
      95          264 :       DO ispin = 1, nspin
      96              :          etr = 0.0_dp
      97          132 :          CALL dbcsr_dot(hmat(1)%matrix, pmat(ispin)%matrix, etr)
      98          264 :          ecore = ecore + etr
      99              :       END DO
     100              : 
     101          132 :       CALL timestop(handle)
     102              : 
     103          132 :    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       311980 :    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       311980 :       CALL timeset(routineN, handle)
     129              : 
     130       311980 :       nc = SIZE(pmat, 2)
     131              : 
     132       311980 :       ecore = 0.0_dp
     133       672542 :       DO ispin = 1, nspin
     134      1775632 :          DO ic = 1, nc
     135              :             etr = 0.0_dp
     136      1103090 :             CALL dbcsr_dot(hmat(1, ic)%matrix, pmat(ispin, ic)%matrix, etr)
     137      1463652 :             ecore = ecore + etr
     138              :          END DO
     139              :       END DO
     140              : 
     141       311980 :       CALL timestop(handle)
     142              : 
     143       311980 :    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           56 :    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              :       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        18551 :    SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, &
     198        18551 :                                       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        18551 :       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        18551 :       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        18551 :       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        18551 :          DIMENSION(:), POINTER                           :: nl_iterator
     222              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     223        18551 :          POINTER                                         :: sab_core
     224        18551 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     225              :       TYPE(qs_energy_type), POINTER                      :: energy
     226        18551 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     227        18551 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     228              :       TYPE(virial_type), POINTER                         :: virial
     229              : 
     230        18551 :       CALL timeset(routineN, handle)
     231              : 
     232        18551 :       NULLIFY (atomic_kind_set)
     233        18551 :       NULLIFY (qs_kind_set)
     234        18551 :       NULLIFY (energy)
     235        18551 :       NULLIFY (atprop)
     236        18551 :       NULLIFY (force)
     237        18551 :       NULLIFY (particle_set)
     238              : 
     239        18551 :       group = para_env
     240              : 
     241        18551 :       only_molecule = .FALSE.
     242        18551 :       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        18551 :                       virial=virial)
     253              : 
     254              :       ! Allocate work storage
     255        18551 :       nkind = SIZE(atomic_kind_set)
     256        18551 :       natom = SIZE(particle_set)
     257              : 
     258        18551 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     259              : 
     260        92755 :       ALLOCATE (alpha(nkind), radius(nkind), zeff(nkind))
     261        51550 :       alpha(:) = 0.0_dp
     262        51550 :       radius(:) = 0.0_dp
     263        51550 :       zeff(:) = 0.0_dp
     264              : 
     265        18551 :       IF (calculate_forces) THEN
     266         6341 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
     267              :       END IF
     268              : 
     269        18551 :       atenergy = .FALSE.
     270        18551 :       IF (ASSOCIATED(atprop)) THEN
     271        18551 :          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        51550 :       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        51550 :                           zeff=zeff(ikind))
     282              :       END DO
     283              : 
     284        18551 :       ecore_overlap = 0.0_dp
     285        18551 :       pv_loc = 0.0_dp
     286              : 
     287        18551 :       CALL neighbor_list_iterator_create(nl_iterator, sab_core)
     288        90111 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     289        71560 :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
     290        71560 :          zab = zeff(ikind)*zeff(jkind)
     291        71560 :          aab = alpha(ikind)*alpha(jkind)/(alpha(ikind) + alpha(jkind))
     292        71560 :          rootaab = SQRT(aab)
     293        71560 :          fab = 2.0_dp*oorootpi*zab*rootaab
     294        71560 :          rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     295        90111 :          IF (rab2 > 1.e-8_dp) THEN
     296        36621 :             IF (ikind == jkind .AND. iatom == jatom) THEN
     297              :                f = 0.5_dp
     298              :             ELSE
     299        36371 :                f = 1.0_dp
     300              :             END IF
     301        36621 :             dab = SQRT(rab2)
     302        36621 :             eab = zab*erfc(rootaab*dab)/dab
     303        36621 :             ecore_overlap = ecore_overlap + f*eab
     304        36621 :             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        36621 :             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        36621 :             IF (calculate_forces) THEN
     313        51292 :                deab(:) = rab(:)*f*(eab + fab*EXP(-aab*rab2))/rab2
     314        12823 :                atom_a = atom_of_kind(iatom)
     315        12823 :                atom_b = atom_of_kind(jatom)
     316        51292 :                force(ikind)%core_overlap(:, atom_a) = force(ikind)%core_overlap(:, atom_a) + deab(:)
     317        51292 :                force(jkind)%core_overlap(:, atom_b) = force(jkind)%core_overlap(:, atom_b) - deab(:)
     318        12823 :                IF (use_virial) THEN
     319          969 :                   CALL virial_pair_force(pv_loc, 1._dp, deab, rab)
     320              :                END IF
     321              :             END IF
     322              :          END IF
     323              :       END DO
     324        18551 :       CALL neighbor_list_iterator_release(nl_iterator)
     325              : 
     326        18551 :       DEALLOCATE (alpha, radius, zeff)
     327        18551 :       IF (calculate_forces) THEN
     328         6341 :          DEALLOCATE (atom_of_kind)
     329              :       END IF
     330        18551 :       IF (calculate_forces .AND. use_virial) THEN
     331         8788 :          virial%pv_ecore_overlap = virial%pv_ecore_overlap + pv_loc
     332         8788 :          virial%pv_virial = virial%pv_virial + pv_loc
     333              :       END IF
     334              : 
     335        18551 :       CALL group%sum(ecore_overlap)
     336              : 
     337        18551 :       energy%core_overlap = ecore_overlap
     338              : 
     339        18551 :       IF (PRESENT(E_overlap_core)) THEN
     340         2066 :          E_overlap_core = energy%core_overlap
     341              :       END IF
     342              : 
     343        18551 :       CALL timestop(handle)
     344              : 
     345        37102 :    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        18119 :    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        18119 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     367              :       TYPE(atprop_type), POINTER                         :: atprop
     368              :       TYPE(distribution_1d_type), POINTER                :: local_particles
     369        18119 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     370              :       TYPE(qs_energy_type), POINTER                      :: energy
     371        18119 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     372              : 
     373              : ! -------------------------------------------------------------------------
     374              : 
     375        18119 :       NULLIFY (atprop)
     376        18119 :       CALL timeset(routineN, handle)
     377              : 
     378              :       CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
     379        18119 :                       qs_kind_set=qs_kind_set, energy=energy, atprop=atprop)
     380              : 
     381        18119 :       ecore_self = 0.0_dp
     382              : 
     383        50636 :       DO ikind = 1, SIZE(atomic_kind_set)
     384        32517 :          CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
     385        32517 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge)
     386        50636 :          ecore_self = ecore_self - REAL(natom, dp)*zeff**2*SQRT(alpha_core_charge)
     387              :       END DO
     388              : 
     389        18119 :       energy%core_self = ecore_self/SQRT(twopi)
     390        18119 :       IF (PRESENT(E_self_core)) THEN
     391         1634 :          E_self_core = energy%core_self
     392              :       END IF
     393              : 
     394        18119 :       IF (ASSOCIATED(atprop)) THEN
     395        18119 :          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        18119 :       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        18119 :       CALL timestop(handle)
     430              : 
     431        18119 :    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
        

Generated by: LCOV version 2.0-1