LCOV - code coverage report
Current view: top level - src - se_fock_matrix_coulomb.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 487 710 68.6 %
Date: 2024-03-28 07:31:50 Functions: 3 4 75.0 %

          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 Module that collects all Coulomb parts of the fock matrix construction
      10             : !>
      11             : !> \author Teodoro Laino (05.2009) [tlaino] - split and module reorganization
      12             : !> \par History
      13             : !>      Teodoro Laino (04.2008) [tlaino] - University of Zurich : d-orbitals
      14             : !>      Teodoro Laino (09.2008) [tlaino] - University of Zurich : Speed-up
      15             : !>      Teodoro Laino (09.2008) [tlaino] - University of Zurich : Periodic SE
      16             : !>      Teodoro Laino (05.2009) [tlaino] - Stress Tensor
      17             : !>
      18             : ! **************************************************************************************************
      19             : MODULE se_fock_matrix_coulomb
      20             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      21             :                                               get_atomic_kind_set
      22             :    USE atprop_types,                    ONLY: atprop_type
      23             :    USE cell_types,                      ONLY: cell_type
      24             :    USE cp_control_types,                ONLY: dft_control_type,&
      25             :                                               semi_empirical_control_type
      26             :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      27             :                                               dbcsr_deallocate_matrix_set
      28             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29             :                                               cp_logger_type
      30             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      31             :                                               cp_print_key_unit_nr
      32             :    USE dbcsr_api,                       ONLY: dbcsr_add,&
      33             :                                               dbcsr_distribute,&
      34             :                                               dbcsr_get_block_diag,&
      35             :                                               dbcsr_get_block_p,&
      36             :                                               dbcsr_p_type,&
      37             :                                               dbcsr_replicate_all,&
      38             :                                               dbcsr_set,&
      39             :                                               dbcsr_sum_replicated
      40             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      41             :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      42             :                                               ewald_environment_type
      43             :    USE ewald_pw_types,                  ONLY: ewald_pw_get,&
      44             :                                               ewald_pw_type
      45             :    USE ewalds_multipole,                ONLY: ewald_multipole_evaluate
      46             :    USE fist_neighbor_list_control,      ONLY: list_control
      47             :    USE fist_nonbond_env_types,          ONLY: fist_nonbond_env_type
      48             :    USE input_constants,                 ONLY: &
      49             :         do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
      50             :         do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1, do_se_IS_slater
      51             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      52             :                                               section_vals_type
      53             :    USE kinds,                           ONLY: dp
      54             :    USE message_passing,                 ONLY: mp_para_env_type
      55             :    USE multipole_types,                 ONLY: do_multipole_charge,&
      56             :                                               do_multipole_dipole,&
      57             :                                               do_multipole_none,&
      58             :                                               do_multipole_quadrupole
      59             :    USE particle_types,                  ONLY: particle_type
      60             :    USE pw_poisson_types,                ONLY: do_ewald_ewald
      61             :    USE qs_energy_types,                 ONLY: qs_energy_type
      62             :    USE qs_environment_types,            ONLY: get_qs_env,&
      63             :                                               qs_environment_type
      64             :    USE qs_force_types,                  ONLY: qs_force_type
      65             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      66             :                                               qs_kind_type
      67             :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      68             :                                               neighbor_list_iterate,&
      69             :                                               neighbor_list_iterator_create,&
      70             :                                               neighbor_list_iterator_p_type,&
      71             :                                               neighbor_list_iterator_release,&
      72             :                                               neighbor_list_set_p_type
      73             :    USE se_fock_matrix_integrals,        ONLY: &
      74             :         dfock2C, dfock2C_r3, dfock2_1el, dfock2_1el_r3, fock2C, fock2C_ew, fock2C_r3, fock2_1el, &
      75             :         fock2_1el_ew, fock2_1el_r3, se_coulomb_ij_interaction
      76             :    USE semi_empirical_int_arrays,       ONLY: rij_threshold,&
      77             :                                               se_orbital_pointer
      78             :    USE semi_empirical_integrals,        ONLY: corecore_el,&
      79             :                                               dcorecore_el
      80             :    USE semi_empirical_mpole_methods,    ONLY: quadrupole_sph_to_cart
      81             :    USE semi_empirical_mpole_types,      ONLY: nddo_mpole_type,&
      82             :                                               semi_empirical_mpole_type
      83             :    USE semi_empirical_store_int_types,  ONLY: semi_empirical_si_type
      84             :    USE semi_empirical_types,            ONLY: get_se_param,&
      85             :                                               se_int_control_type,&
      86             :                                               se_taper_type,&
      87             :                                               semi_empirical_p_type,&
      88             :                                               semi_empirical_type,&
      89             :                                               setup_se_int_control_type
      90             :    USE semi_empirical_utils,            ONLY: finalize_se_taper,&
      91             :                                               get_se_type,&
      92             :                                               initialize_se_taper
      93             :    USE virial_methods,                  ONLY: virial_pair_force
      94             :    USE virial_types,                    ONLY: virial_type
      95             : #include "./base/base_uses.f90"
      96             : 
      97             :    IMPLICIT NONE
      98             :    PRIVATE
      99             : 
     100             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_coulomb'
     101             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
     102             : 
     103             :    PUBLIC :: build_fock_matrix_coulomb, build_fock_matrix_coulomb_lr, &
     104             :              build_fock_matrix_coul_lr_r3
     105             : 
     106             : CONTAINS
     107             : 
     108             : ! **************************************************************************************************
     109             : !> \brief Construction of the Coulomb part of the Fock matrix
     110             : !> \param qs_env ...
     111             : !> \param ks_matrix ...
     112             : !> \param matrix_p ...
     113             : !> \param energy ...
     114             : !> \param calculate_forces ...
     115             : !> \param store_int_env ...
     116             : !> \author JGH
     117             : ! **************************************************************************************************
     118       41070 :    SUBROUTINE build_fock_matrix_coulomb(qs_env, ks_matrix, matrix_p, energy, calculate_forces, &
     119             :                                         store_int_env)
     120             : 
     121             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     122             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix, matrix_p
     123             :       TYPE(qs_energy_type), POINTER                      :: energy
     124             :       LOGICAL, INTENT(in)                                :: calculate_forces
     125             :       TYPE(semi_empirical_si_type), POINTER              :: store_int_env
     126             : 
     127             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coulomb'
     128             : 
     129             :       INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, ispin, itype, jatom, jkind, &
     130             :          natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nspins
     131       41070 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, se_natorb
     132             :       LOGICAL                                            :: anag, atener, check, defined, found, &
     133             :                                                             switch, use_virial
     134       41070 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: se_defined
     135             :       REAL(KIND=dp)                                      :: delta, dr1, ecore2, ecores
     136             :       REAL(KIND=dp), DIMENSION(2)                        :: ecab
     137             :       REAL(KIND=dp), DIMENSION(2025)                     :: pa_a, pa_b, pb_a, pb_b
     138             :       REAL(KIND=dp), DIMENSION(3)                        :: force_ab, rij
     139             :       REAL(KIND=dp), DIMENSION(45, 45)                   :: p_block_tot_a, p_block_tot_b
     140       41070 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ksa_block_a, ksa_block_b, ksb_block_a, &
     141       41070 :                                                             ksb_block_b, pa_block_a, pa_block_b, &
     142       41070 :                                                             pb_block_a, pb_block_b
     143       41070 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     144             :       TYPE(atprop_type), POINTER                         :: atprop
     145             :       TYPE(cell_type), POINTER                           :: cell
     146       41070 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: diagmat_ks, diagmat_p
     147             :       TYPE(dft_control_type), POINTER                    :: dft_control
     148             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     149             :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     150             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     151             :       TYPE(neighbor_list_iterator_p_type), &
     152       41070 :          DIMENSION(:), POINTER                           :: nl_iterator
     153             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     154       41070 :          POINTER                                         :: sab_se
     155       41070 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     156       41070 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     157       41070 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     158             :       TYPE(se_int_control_type)                          :: se_int_control
     159             :       TYPE(se_taper_type), POINTER                       :: se_taper
     160             :       TYPE(semi_empirical_control_type), POINTER         :: se_control
     161       41070 :       TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
     162             :       TYPE(semi_empirical_type), POINTER                 :: se_kind_a, se_kind_b
     163             :       TYPE(virial_type), POINTER                         :: virial
     164             : 
     165       41070 :       CALL timeset(routineN, handle)
     166             : 
     167       41070 :       NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, &
     168       41070 :                se_control, se_taper, virial, atprop)
     169             : 
     170             :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, &
     171             :                       para_env=para_env, sab_se=sab_se, atomic_kind_set=atomic_kind_set, atprop=atprop, &
     172       41070 :                       qs_kind_set=qs_kind_set, particle_set=particle_set, virial=virial)
     173             : 
     174             :       ! Parameters
     175       41070 :       CALL initialize_se_taper(se_taper, coulomb=.TRUE.)
     176       41070 :       se_control => dft_control%qs_control%se_control
     177       41070 :       anag = se_control%analytical_gradients
     178       41070 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     179       41070 :       atener = atprop%energy
     180             : 
     181             :       CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, &
     182             :                                      do_ewald_gks=se_control%do_ewald_gks, integral_screening=se_control%integral_screening, &
     183             :                                      shortrange=(se_control%do_ewald .OR. se_control%do_ewald_gks), &
     184       80404 :                                      max_multipole=se_control%max_multipole, pc_coulomb_int=.FALSE.)
     185       41070 :       IF (se_control%do_ewald_gks) THEN
     186         106 :          CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
     187         106 :          CALL ewald_env_get(ewald_env, alpha=se_int_control%ewald_gks%alpha)
     188             :          CALL ewald_pw_get(ewald_pw, pw_big_pool=se_int_control%ewald_gks%pw_pool, &
     189         106 :                            dg=se_int_control%ewald_gks%dg)
     190             :       END IF
     191             : 
     192       41070 :       nspins = dft_control%nspins
     193       41070 :       CPASSERT(ASSOCIATED(matrix_p))
     194       41070 :       CPASSERT(SIZE(ks_matrix) > 0)
     195             : 
     196       41070 :       nkind = SIZE(atomic_kind_set)
     197       41070 :       IF (calculate_forces) THEN
     198        3002 :          CALL get_qs_env(qs_env=qs_env, force=force)
     199        3002 :          delta = se_control%delta
     200        3002 :          CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     201             :       END IF
     202             : 
     203       41070 :       CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
     204       41070 :       CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
     205             : 
     206       84704 :       DO ispin = 1, nspins
     207             :          ! Allocate diagonal block matrices
     208       43634 :          ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
     209       43634 :          CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
     210       43634 :          CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
     211       43634 :          CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
     212       43634 :          CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
     213       84704 :          CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
     214             :       END DO
     215             : 
     216       41070 :       ecore2 = 0.0_dp
     217       41070 :       itype = get_se_type(dft_control%qs_control%method_id)
     218      287490 :       ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
     219      143288 :       DO ikind = 1, nkind
     220      102218 :          CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
     221      102218 :          se_kind_list(ikind)%se_param => se_kind_a
     222      102218 :          CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
     223      102218 :          se_defined(ikind) = (defined .AND. natorb_a >= 1)
     224      245506 :          se_natorb(ikind) = natorb_a
     225             :       END DO
     226             : 
     227       41070 :       CALL neighbor_list_iterator_create(nl_iterator, sab_se)
     228     9947992 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     229     9906922 :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
     230     9906922 :          IF (.NOT. se_defined(ikind)) CYCLE
     231     9906922 :          IF (.NOT. se_defined(jkind)) CYCLE
     232     9906922 :          se_kind_a => se_kind_list(ikind)%se_param
     233     9906922 :          se_kind_b => se_kind_list(jkind)%se_param
     234     9906922 :          natorb_a = se_natorb(ikind)
     235     9906922 :          natorb_b = se_natorb(jkind)
     236     9906922 :          natorb_a2 = natorb_a**2
     237     9906922 :          natorb_b2 = natorb_b**2
     238             : 
     239     9906922 :          IF (inode == 1) THEN
     240             :             CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
     241      639804 :                                    row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
     242      639804 :             CPASSERT(ASSOCIATED(pa_block_a))
     243      639804 :             check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
     244           0 :             CPASSERT(check)
     245             :             CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
     246      639804 :                                    row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
     247      639804 :             CPASSERT(ASSOCIATED(ksa_block_a))
     248    10124040 :             p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
     249     1279608 :             pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/))
     250      639804 :             IF (nspins == 2) THEN
     251             :                CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
     252        8389 :                                       row=iatom, col=iatom, BLOCK=pa_block_b, found=found)
     253        8389 :                CPASSERT(ASSOCIATED(pa_block_b))
     254        8389 :                check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
     255           0 :                CPASSERT(check)
     256             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
     257        8389 :                                       row=iatom, col=iatom, BLOCK=ksa_block_b, found=found)
     258        8389 :                CPASSERT(ASSOCIATED(ksa_block_b))
     259      117595 :                p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
     260       16778 :                pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/))
     261             :             END IF
     262             : 
     263             :          END IF
     264             : 
     265    39627688 :          dr1 = DOT_PRODUCT(rij, rij)
     266     9947992 :          IF (dr1 > rij_threshold) THEN
     267             :             ! Determine the order of the atoms, and in case switch them..
     268     9710345 :             IF (iatom <= jatom) THEN
     269     5102392 :                switch = .FALSE.
     270             :             ELSE
     271     4607953 :                switch = .TRUE.
     272             :             END IF
     273             :             ! Retrieve blocks for KS and P
     274             :             CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
     275     9710345 :                                    row=jatom, col=jatom, BLOCK=pb_block_a, found=found)
     276     9710345 :             CPASSERT(ASSOCIATED(pb_block_a))
     277     9710345 :             check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
     278           0 :             CPASSERT(check)
     279             :             CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
     280     9710345 :                                    row=jatom, col=jatom, BLOCK=ksb_block_a, found=found)
     281     9710345 :             CPASSERT(ASSOCIATED(ksb_block_a))
     282   148337719 :             p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
     283    19420690 :             pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/))
     284             :             ! Handle more than one configuration
     285     9710345 :             IF (nspins == 2) THEN
     286             :                CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
     287      225262 :                                       row=jatom, col=jatom, BLOCK=pb_block_b, found=found)
     288      225262 :                CPASSERT(ASSOCIATED(pb_block_b))
     289      225262 :                check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
     290           0 :                CPASSERT(check)
     291             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
     292      225262 :                                       row=jatom, col=jatom, BLOCK=ksb_block_b, found=found)
     293      225262 :                CPASSERT(ASSOCIATED(ksb_block_b))
     294      225262 :                check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
     295           0 :                CPASSERT(check)
     296     2764564 :                p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
     297      450524 :                pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/))
     298             :             END IF
     299             : 
     300    19420690 :             SELECT CASE (dft_control%qs_control%method_id)
     301             :             CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
     302             :                   do_method_rm1, do_method_mndod, do_method_pnnl)
     303             : 
     304             :                ! Two-centers One-electron terms
     305     9710345 :                IF (nspins == 1) THEN
     306     9485083 :                   ecab = 0._dp
     307             :                   CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
     308             :                                  pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
     309     9485083 :                                  se_taper=se_taper, store_int_env=store_int_env)
     310     9485083 :                   ecore2 = ecore2 + ecab(1) + ecab(2)
     311      225262 :                ELSE IF (nspins == 2) THEN
     312      225262 :                   ecab = 0._dp
     313             :                   CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
     314             :                                  pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag, &
     315      225262 :                                  se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
     316             :                   CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_b, ksb_block_b, &
     317             :                                  pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
     318      225262 :                                  se_taper=se_taper, store_int_env=store_int_env)
     319      225262 :                   ecore2 = ecore2 + ecab(1) + ecab(2)
     320             :                END IF
     321     9710345 :                IF (atener) THEN
     322         765 :                   atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1)
     323         765 :                   atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2)
     324             :                END IF
     325             :                ! Coulomb Terms
     326     9710345 :                IF (nspins == 1) THEN
     327             :                   CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
     328             :                               ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
     329     9485083 :                               store_int_env=store_int_env)
     330      225262 :                ELSE IF (nspins == 2) THEN
     331             :                   CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
     332             :                               ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
     333      225262 :                               store_int_env=store_int_env)
     334             : 
     335             :                   CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, &
     336             :                               ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
     337      225262 :                               store_int_env=store_int_env)
     338             :                END IF
     339             : 
     340     9710345 :                IF (calculate_forces) THEN
     341      302154 :                   atom_a = atom_of_kind(iatom)
     342      302154 :                   atom_b = atom_of_kind(jatom)
     343             : 
     344             :                   ! Derivatives of the Two-centre One-electron terms
     345      302154 :                   force_ab = 0.0_dp
     346      302154 :                   IF (nspins == 1) THEN
     347             :                      CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_a, pb_a, itype=itype, anag=anag, &
     348             :                                      se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
     349      301993 :                                      delta=delta)
     350         161 :                   ELSE IF (nspins == 2) THEN
     351             :                      CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_block_a, pb_block_a, itype=itype, anag=anag, &
     352         161 :                                      se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
     353             :                      CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_b, pb_b, itype=itype, anag=anag, &
     354         161 :                                      se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
     355             :                   END IF
     356      302154 :                   IF (use_virial) THEN
     357       18728 :                      CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
     358             :                   END IF
     359             : 
     360             :                   ! Sum up force components
     361      302154 :                   force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
     362      302154 :                   force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
     363             : 
     364      302154 :                   force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
     365      302154 :                   force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
     366             : 
     367      302154 :                   force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
     368      302154 :                   force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
     369             : 
     370             :                   ! Derivatives of the Coulomb Terms
     371      302154 :                   force_ab = 0._dp
     372      302154 :                   IF (nspins == 1) THEN
     373             :                      CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, &
     374      301993 :                                   anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
     375         161 :                   ELSE IF (nspins == 2) THEN
     376             :                      CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
     377         161 :                                   anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
     378             : 
     379             :                      CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
     380         161 :                                   anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
     381             :                   END IF
     382      302154 :                   IF (switch) THEN
     383      149436 :                      force_ab(1) = -force_ab(1)
     384      149436 :                      force_ab(2) = -force_ab(2)
     385      149436 :                      force_ab(3) = -force_ab(3)
     386             :                   END IF
     387      302154 :                   IF (use_virial) THEN
     388       18728 :                      CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
     389             :                   END IF
     390             :                   ! Sum up force components
     391      302154 :                   force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
     392      302154 :                   force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
     393             : 
     394      302154 :                   force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
     395      302154 :                   force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
     396             : 
     397      302154 :                   force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
     398      302154 :                   force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
     399             :                END IF
     400             :             CASE DEFAULT
     401     9710345 :                CPABORT("")
     402             :             END SELECT
     403             :          ELSE
     404      196577 :             IF (se_int_control%do_ewald_gks) THEN
     405         159 :                CPASSERT(iatom == jatom)
     406             :                ! Two-centers One-electron terms
     407         159 :                ecores = 0._dp
     408         159 :                IF (nspins == 1) THEN
     409             :                   CALL fock2_1el_ew(se_kind_a, rij, ksa_block_a, pa_a, &
     410             :                                     ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, &
     411         159 :                                     se_taper=se_taper, store_int_env=store_int_env)
     412           0 :                ELSE IF (nspins == 2) THEN
     413             :                   CALL fock2_1el_ew(se_kind_a, rij, ksa_block_a, pa_block_a, &
     414             :                                     ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, &
     415           0 :                                     se_taper=se_taper, store_int_env=store_int_env)
     416             :                   CALL fock2_1el_ew(se_kind_a, rij, ksa_block_b, pa_b, &
     417             :                                     ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, &
     418           0 :                                     se_taper=se_taper, store_int_env=store_int_env)
     419             :                END IF
     420         159 :                ecore2 = ecore2 + ecores
     421         159 :                IF (atener) THEN
     422           0 :                   atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecores
     423             :                END IF
     424             :                ! Coulomb Terms
     425         159 :                IF (nspins == 1) THEN
     426             :                   CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a, &
     427             :                                  factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
     428         159 :                                  store_int_env=store_int_env)
     429           0 :                ELSE IF (nspins == 2) THEN
     430             :                   CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a, &
     431             :                                  factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
     432           0 :                                  store_int_env=store_int_env)
     433             :                   CALL fock2C_ew(se_kind_a, rij, p_block_tot_a, ksa_block_b, &
     434             :                                  factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
     435           0 :                                  store_int_env=store_int_env)
     436             :                END IF
     437             :             END IF
     438             :          END IF
     439             :       END DO
     440       41070 :       CALL neighbor_list_iterator_release(nl_iterator)
     441             : 
     442       41070 :       DEALLOCATE (se_kind_list, se_defined, se_natorb)
     443             : 
     444       84704 :       DO ispin = 1, nspins
     445       43634 :          CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
     446       43634 :          CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
     447             :          CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
     448       84704 :                         1.0_dp, 1.0_dp)
     449             :       END DO
     450       41070 :       CALL dbcsr_deallocate_matrix_set(diagmat_p)
     451       41070 :       CALL dbcsr_deallocate_matrix_set(diagmat_ks)
     452             : 
     453             :       ! Two-centers one-electron terms
     454       41070 :       CALL para_env%sum(ecore2)
     455       41070 :       energy%hartree = ecore2 - energy%core
     456             : !      WRITE ( *, * ) 'IN SE_F_COUL', ecore2, energy%core
     457             : 
     458       41070 :       CALL finalize_se_taper(se_taper)
     459             : 
     460       41070 :       CALL timestop(handle)
     461       82140 :    END SUBROUTINE build_fock_matrix_coulomb
     462             : 
     463             : ! **************************************************************************************************
     464             : !> \brief  Long-Range part for SE Coulomb interactions
     465             : !> \param qs_env ...
     466             : !> \param ks_matrix ...
     467             : !> \param matrix_p ...
     468             : !> \param energy ...
     469             : !> \param calculate_forces ...
     470             : !> \param store_int_env ...
     471             : !> \date   08.2008 [created]
     472             : !> \author Teodoro Laino [tlaino] - University of Zurich
     473             : ! **************************************************************************************************
     474        1630 :    SUBROUTINE build_fock_matrix_coulomb_lr(qs_env, ks_matrix, matrix_p, energy, &
     475             :                                            calculate_forces, store_int_env)
     476             : 
     477             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     478             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix, matrix_p
     479             :       TYPE(qs_energy_type), POINTER                      :: energy
     480             :       LOGICAL, INTENT(in)                                :: calculate_forces
     481             :       TYPE(semi_empirical_si_type), POINTER              :: store_int_env
     482             : 
     483             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coulomb_lr'
     484             : 
     485             :       INTEGER :: atom_a, ewald_type, forces_g_size, handle, iatom, ikind, ilist, indi, indj, &
     486             :          ispin, itype, iw, jint, natoms, natorb_a, nkind, nlocal_particles, node, nparticle_local, &
     487             :          nspins, size_1c_int
     488        1630 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     489             :       LOGICAL                                            :: anag, atener, defined, found, use_virial
     490             :       LOGICAL, DIMENSION(3)                              :: task
     491             :       REAL(KIND=dp)                                      :: e_neut, e_self, energy_glob, &
     492             :                                                             energy_local, enuc, fac, tmp
     493        1630 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: forces_g, forces_r
     494             :       REAL(KIND=dp), DIMENSION(3)                        :: force_a
     495             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_glob, pv_local, qcart
     496             :       REAL(KIND=dp), DIMENSION(5)                        :: qsph
     497        1630 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ksa_block_a, pa_block_a
     498        1630 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     499             :       TYPE(atprop_type), POINTER                         :: atprop
     500             :       TYPE(cell_type), POINTER                           :: cell
     501             :       TYPE(cp_logger_type), POINTER                      :: logger
     502        1630 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: diagmat_ks, diagmat_p
     503             :       TYPE(dft_control_type), POINTER                    :: dft_control
     504             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     505             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     506             :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     507             :       TYPE(fist_nonbond_env_type), POINTER               :: se_nonbond_env
     508             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     509             :       TYPE(nddo_mpole_type), POINTER                     :: se_nddo_mpole
     510        1630 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     511        1630 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     512        1630 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     513             :       TYPE(section_vals_type), POINTER                   :: se_section
     514             :       TYPE(semi_empirical_control_type), POINTER         :: se_control
     515             :       TYPE(semi_empirical_mpole_type), POINTER           :: mpole
     516             :       TYPE(semi_empirical_type), POINTER                 :: se_kind_a
     517             :       TYPE(virial_type), POINTER                         :: virial
     518             : 
     519        1630 :       CALL timeset(routineN, handle)
     520             : 
     521        1630 :       NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, local_particles, &
     522        1630 :                se_control, ewald_env, ewald_pw, se_nddo_mpole, se_nonbond_env, se_section, mpole, &
     523        1630 :                logger, virial, atprop)
     524             : 
     525        1630 :       logger => cp_get_default_logger()
     526             :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env, &
     527             :                       atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set, &
     528             :                       ewald_env=ewald_env, &
     529             :                       local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole, &
     530        1630 :                       se_nonbond_env=se_nonbond_env, virial=virial, atprop=atprop)
     531             : 
     532        4538 :       nlocal_particles = SUM(local_particles%n_el(:))
     533        1630 :       natoms = SIZE(particle_set)
     534        1630 :       CALL ewald_env_get(ewald_env, ewald_type=ewald_type)
     535             :       SELECT CASE (ewald_type)
     536             :       CASE (do_ewald_ewald)
     537           0 :          forces_g_size = nlocal_particles
     538             :       CASE DEFAULT
     539        1630 :          CPABORT("Periodic SE implemented only for standard EWALD sums.")
     540             :       END SELECT
     541             : 
     542             :       ! Parameters
     543        1630 :       se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE")
     544        1630 :       se_control => dft_control%qs_control%se_control
     545        1630 :       anag = se_control%analytical_gradients
     546        1630 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
     547        1630 :       atener = atprop%energy
     548             : 
     549        1630 :       nspins = dft_control%nspins
     550        1630 :       CPASSERT(ASSOCIATED(matrix_p))
     551        1630 :       CPASSERT(SIZE(ks_matrix) > 0)
     552             : 
     553        1630 :       CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
     554        1630 :       CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
     555             : 
     556        1630 :       nkind = SIZE(atomic_kind_set)
     557        3260 :       DO ispin = 1, nspins
     558             :          ! Allocate diagonal block matrices
     559        1630 :          ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
     560        1630 :          CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
     561        1630 :          CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
     562        1630 :          CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
     563        1630 :          CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
     564        3260 :          CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
     565             :       END DO
     566             : 
     567             :       ! Check for implemented SE methods
     568        1630 :       SELECT CASE (dft_control%qs_control%method_id)
     569             :       CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
     570             :             do_method_rm1, do_method_mndod, do_method_pnnl)
     571           0 :          itype = get_se_type(dft_control%qs_control%method_id)
     572             :       CASE DEFAULT
     573        1630 :          CPABORT("")
     574             :       END SELECT
     575             : 
     576             :       ! Zero arrays and possibly build neighbor lists
     577        1630 :       energy_local = 0.0_dp
     578        1630 :       energy_glob = 0.0_dp
     579        1630 :       e_neut = 0.0_dp
     580        1630 :       e_self = 0.0_dp
     581        1630 :       task = .FALSE.
     582        1630 :       SELECT CASE (se_control%max_multipole)
     583             :       CASE (do_multipole_none)
     584             :          ! Do Nothing
     585             :       CASE (do_multipole_charge)
     586           0 :          task(1) = .TRUE.
     587             :       CASE (do_multipole_dipole)
     588           0 :          task = .TRUE.
     589           0 :          task(3) = .FALSE.
     590             :       CASE (do_multipole_quadrupole)
     591        4832 :          task = .TRUE.
     592             :       CASE DEFAULT
     593        1630 :          CPABORT("")
     594             :       END SELECT
     595             : 
     596             :       ! Build-up neighbor lists for real-space part of Ewald multipoles
     597             :       CALL list_control(atomic_kind_set, particle_set, local_particles, &
     598        1630 :                         cell, se_nonbond_env, para_env, se_section)
     599             : 
     600        1630 :       enuc = 0.0_dp
     601        1630 :       energy%core_overlap = 0.0_dp
     602       16610 :       se_nddo_mpole%charge = 0.0_dp
     603       61550 :       se_nddo_mpole%dipole = 0.0_dp
     604      196370 :       se_nddo_mpole%quadrupole = 0.0_dp
     605             : 
     606        3260 :       DO ispin = 1, nspins
     607             :          ! Compute the NDDO mpole expansion
     608        6168 :          DO ikind = 1, nkind
     609        2908 :             CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
     610        2908 :             CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
     611             : 
     612        2908 :             IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     613             : 
     614        2908 :             nparticle_local = local_particles%n_el(ikind)
     615       14936 :             DO ilist = 1, nparticle_local
     616        7490 :                iatom = local_particles%list(ikind)%array(ilist)
     617             :                CALL dbcsr_get_block_p(matrix=diagmat_p(ispin)%matrix, &
     618        7490 :                                       row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
     619        7490 :                CPASSERT(ASSOCIATED(pa_block_a))
     620             :                ! Nuclei
     621        7490 :                IF (task(1) .AND. ispin == 1) se_nddo_mpole%charge(iatom) = se_kind_a%zeff
     622             :                ! Electrons
     623        7490 :                size_1c_int = SIZE(se_kind_a%w_mpole)
     624       70294 :                DO jint = 1, size_1c_int
     625       62804 :                   mpole => se_kind_a%w_mpole(jint)%mpole
     626       62804 :                   indi = se_orbital_pointer(mpole%indi)
     627       62804 :                   indj = se_orbital_pointer(mpole%indj)
     628       62804 :                   fac = 1.0_dp
     629       62804 :                   IF (indi /= indj) fac = 2.0_dp
     630             : 
     631             :                   ! Charge
     632       62804 :                   IF (mpole%task(1) .AND. task(1)) THEN
     633             :                      se_nddo_mpole%charge(iatom) = se_nddo_mpole%charge(iatom) + &
     634       18092 :                                                    fac*pa_block_a(indi, indj)*mpole%c
     635             :                   END IF
     636             : 
     637             :                   ! Dipole
     638       62804 :                   IF (mpole%task(2) .AND. task(2)) THEN
     639             :                      se_nddo_mpole%dipole(:, iatom) = se_nddo_mpole%dipole(:, iatom) + &
     640       89019 :                                                       fac*pa_block_a(indi, indj)*mpole%d(:)
     641             :                   END IF
     642             : 
     643             :                   ! Quadrupole
     644       70294 :                   IF (mpole%task(3) .AND. task(3)) THEN
     645      151740 :                      qsph = fac*mpole%qs*pa_block_a(indi, indj)
     646       25290 :                      CALL quadrupole_sph_to_cart(qcart, qsph)
     647             :                      se_nddo_mpole%quadrupole(:, :, iatom) = se_nddo_mpole%quadrupole(:, :, iatom) + &
     648      328770 :                                                              qcart
     649             :                   END IF
     650             :                END DO
     651             :                ! Print some info about charge, dipole and quadrupole (debug purpose only)
     652       10398 :                IF (debug_this_module) THEN
     653             :                   WRITE (*, '(I5,F12.6,5X,3F12.6,5X,9F12.6)') iatom, se_nddo_mpole%charge(iatom), &
     654             :                      se_nddo_mpole%dipole(:, iatom), se_nddo_mpole%quadrupole(:, :, iatom)
     655             :                END IF
     656             :             END DO
     657             :          END DO
     658             :       END DO
     659       31590 :       CALL para_env%sum(se_nddo_mpole%charge)
     660      121470 :       CALL para_env%sum(se_nddo_mpole%dipole)
     661      391110 :       CALL para_env%sum(se_nddo_mpole%quadrupole)
     662             : 
     663             :       ! Initialize for virial
     664        1630 :       IF (use_virial) THEN
     665           2 :          pv_glob = 0.0_dp
     666           2 :          pv_local = 0.0_dp
     667             :       END IF
     668             : 
     669             :       ! Ewald Multipoles Sum
     670        1630 :       iw = cp_print_key_unit_nr(logger, se_section, "PRINT%EWALD_INFO", extension=".seLog")
     671        1630 :       IF (calculate_forces) THEN
     672          20 :          CALL get_qs_env(qs_env=qs_env, force=force)
     673          20 :          CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     674             : 
     675             :          ! Allocate and zeroing arrays
     676          59 :          ALLOCATE (forces_g(3, forces_g_size))
     677          60 :          ALLOCATE (forces_r(3, natoms))
     678         668 :          forces_g = 0.0_dp
     679        1316 :          forces_r = 0.0_dp
     680             :          CALL ewald_multipole_evaluate( &
     681             :             ewald_env, ewald_pw, se_nonbond_env, cell, &
     682             :             particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, &
     683             :             do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=use_virial, do_efield=.TRUE., &
     684             :             charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, &
     685             :             forces_local=forces_g, forces_glob=forces_r, pv_glob=pv_glob, pv_local=pv_local, &
     686             :             efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, iw=iw, &
     687          20 :             do_debug=.TRUE.)
     688             :          ! Only SR force have to be summed up.. the one in g-space are already fully local..
     689          20 :          CALL para_env%sum(forces_r)
     690             :       ELSE
     691             :          CALL ewald_multipole_evaluate( &
     692             :             ewald_env, ewald_pw, se_nonbond_env, cell, &
     693             :             particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, &
     694             :             do_correction_bonded=.FALSE., do_forces=.FALSE., do_stress=.FALSE., do_efield=.TRUE., &
     695             :             charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, &
     696             :             efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, &
     697        1610 :             iw=iw, do_debug=.TRUE.)
     698             :       END IF
     699        1630 :       CALL cp_print_key_finished_output(iw, logger, se_section, "PRINT%EWALD_INFO")
     700             : 
     701             :       ! Apply correction only when the Integral Scheme is different from Slater
     702        1630 :       IF ((se_control%integral_screening /= do_se_IS_slater) .AND. (.NOT. debug_this_module)) THEN
     703             :          CALL build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, calculate_forces, &
     704             :                                          store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, virial, &
     705        1400 :                                          pv_glob)
     706             :       END IF
     707             : 
     708             :       ! Virial for the long-range part and correction
     709        1630 :       IF (use_virial) THEN
     710             :          ! Sum up contribution of pv_glob on each thread and keep only one copy of pv_local
     711          26 :          virial%pv_virial = virial%pv_virial + pv_glob
     712           2 :          IF (para_env%is_source()) THEN
     713          13 :             virial%pv_virial = virial%pv_virial + pv_local
     714             :          END IF
     715             :       END IF
     716             : 
     717             :       ! Debug Statements
     718             :       IF (debug_this_module) THEN
     719             :          CALL para_env%sum(energy_glob)
     720             :          WRITE (*, *) "TOTAL ENERGY AFTER EWALD:", energy_local + energy_glob + e_neut + e_self, &
     721             :             energy_local, energy_glob, e_neut, e_self
     722             :       END IF
     723             : 
     724             :       ! Modify the KS matrix and possibly compute derivatives
     725             :       node = 0
     726        4538 :       DO ikind = 1, nkind
     727        2908 :          CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
     728        2908 :          CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
     729             : 
     730        2908 :          IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     731             : 
     732        2908 :          nparticle_local = local_particles%n_el(ikind)
     733       14936 :          DO ilist = 1, nparticle_local
     734        7490 :             node = node + 1
     735        7490 :             iatom = local_particles%list(ikind)%array(ilist)
     736       14980 :             DO ispin = 1, nspins
     737             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(ispin)%matrix, &
     738        7490 :                                       row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
     739        7490 :                CPASSERT(ASSOCIATED(ksa_block_a))
     740             : 
     741             :                ! Modify Hamiltonian Matrix accordingly potential, field and electric field gradient
     742        7490 :                size_1c_int = SIZE(se_kind_a%w_mpole)
     743       85274 :                DO jint = 1, size_1c_int
     744       62804 :                   tmp = 0.0_dp
     745       62804 :                   mpole => se_kind_a%w_mpole(jint)%mpole
     746       62804 :                   indi = se_orbital_pointer(mpole%indi)
     747       62804 :                   indj = se_orbital_pointer(mpole%indj)
     748             : 
     749             :                   ! Charge
     750       62804 :                   IF (mpole%task(1) .AND. task(1)) THEN
     751       18092 :                      tmp = tmp + mpole%c*se_nddo_mpole%efield0(iatom)
     752             :                   END IF
     753             : 
     754             :                   ! Dipole
     755       62804 :                   IF (mpole%task(2) .AND. task(2)) THEN
     756       50868 :                      tmp = tmp - DOT_PRODUCT(mpole%d, se_nddo_mpole%efield1(:, iatom))
     757             :                   END IF
     758             : 
     759             :                   ! Quadrupole
     760       62804 :                   IF (mpole%task(3) .AND. task(3)) THEN
     761      328770 :                      tmp = tmp - (1.0_dp/3.0_dp)*SUM(mpole%qc*RESHAPE(se_nddo_mpole%efield2(:, iatom), (/3, 3/)))
     762             :                   END IF
     763             : 
     764       62804 :                   ksa_block_a(indi, indj) = ksa_block_a(indi, indj) + tmp
     765       70294 :                   ksa_block_a(indj, indi) = ksa_block_a(indi, indj)
     766             :                END DO
     767             :             END DO
     768             : 
     769             :             ! Nuclear term and forces
     770        7490 :             IF (task(1)) enuc = enuc + se_kind_a%zeff*se_nddo_mpole%efield0(iatom)
     771        7490 :             IF (atener) THEN
     772             :                atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
     773           0 :                                        0.5_dp*se_kind_a%zeff*se_nddo_mpole%efield0(iatom)
     774             :             END IF
     775       10398 :             IF (calculate_forces) THEN
     776         162 :                atom_a = atom_of_kind(iatom)
     777         648 :                force_a = forces_r(1:3, iatom) + forces_g(1:3, node)
     778             :                ! Derivatives of the periodic Coulomb Terms
     779         648 :                force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_a(:)
     780             :             END IF
     781             :          END DO
     782             :       END DO
     783             :       ! Sum nuclear energy contribution
     784        1630 :       CALL para_env%sum(enuc)
     785        1630 :       energy%core_overlap = energy%core_overlap + energy%core_overlap0 + 0.5_dp*enuc
     786             : 
     787             :       ! Debug Statements
     788             :       IF (debug_this_module) THEN
     789             :          WRITE (*, *) "ENUC: ", enuc*0.5_dp
     790             :       END IF
     791             : 
     792        3260 :       DO ispin = 1, nspins
     793        1630 :          CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
     794        1630 :          CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
     795             :          CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
     796        3260 :                         1.0_dp, 1.0_dp)
     797             :       END DO
     798        1630 :       CALL dbcsr_deallocate_matrix_set(diagmat_p)
     799        1630 :       CALL dbcsr_deallocate_matrix_set(diagmat_ks)
     800             : 
     801             :       ! Set the Fock matrix contribution to SCP
     802        1630 :       IF (calculate_forces) THEN
     803          20 :          DEALLOCATE (forces_g)
     804          20 :          DEALLOCATE (forces_r)
     805             :       END IF
     806             : 
     807        1630 :       CALL timestop(handle)
     808        3260 :    END SUBROUTINE build_fock_matrix_coulomb_lr
     809             : 
     810             : ! **************************************************************************************************
     811             : !> \brief When doing long-range SE calculation, this module computes the correction
     812             : !>        between the mismatch of point-like multipoles and multipoles represented
     813             : !>        with charges
     814             : !> \param qs_env ...
     815             : !> \param ks_matrix ...
     816             : !> \param matrix_p ...
     817             : !> \param energy ...
     818             : !> \param calculate_forces ...
     819             : !> \param store_int_env ...
     820             : !> \param se_nddo_mpole ...
     821             : !> \param task ...
     822             : !> \param diagmat_p ...
     823             : !> \param diagmat_ks ...
     824             : !> \param virial ...
     825             : !> \param pv_glob ...
     826             : !> \author Teodoro Laino [tlaino] - 05.2009
     827             : ! **************************************************************************************************
     828        1400 :    SUBROUTINE build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, &
     829             :                                          calculate_forces, store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, &
     830             :                                          virial, pv_glob)
     831             : 
     832             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     833             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix, matrix_p
     834             :       TYPE(qs_energy_type), POINTER                      :: energy
     835             :       LOGICAL, INTENT(in)                                :: calculate_forces
     836             :       TYPE(semi_empirical_si_type), POINTER              :: store_int_env
     837             :       TYPE(nddo_mpole_type), POINTER                     :: se_nddo_mpole
     838             :       LOGICAL, DIMENSION(3), INTENT(IN)                  :: task
     839             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: diagmat_p, diagmat_ks
     840             :       TYPE(virial_type), POINTER                         :: virial
     841             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: pv_glob
     842             : 
     843             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coul_lrc'
     844             : 
     845             :       INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, itype, jatom, jkind, natorb_a, &
     846             :          natorb_a2, natorb_b, natorb_b2, nkind, nspins, size1, size2
     847        1400 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, se_natorb
     848             :       LOGICAL                                            :: anag, atener, check, defined, found, &
     849             :                                                             switch, use_virial
     850        1400 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: se_defined
     851             :       REAL(KIND=dp)                                      :: delta, dr1, ecore2, enuc, enuclear, &
     852             :                                                             ptens11, ptens12, ptens13, ptens21, &
     853             :                                                             ptens22, ptens23, ptens31, ptens32, &
     854             :                                                             ptens33
     855             :       REAL(KIND=dp), DIMENSION(2)                        :: ecab
     856             :       REAL(KIND=dp), DIMENSION(2025)                     :: pa_a, pa_b, pb_a, pb_b
     857             :       REAL(KIND=dp), DIMENSION(3)                        :: force_ab, force_ab0, rij
     858             :       REAL(KIND=dp), DIMENSION(45, 45)                   :: p_block_tot_a, p_block_tot_b
     859        1400 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: efield0
     860        1400 :       REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2, ksa_block_a, ksa_block_b, &
     861        1400 :          ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, pb_block_a, pb_block_b
     862        1400 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     863             :       TYPE(atprop_type), POINTER                         :: atprop
     864             :       TYPE(cell_type), POINTER                           :: cell
     865             :       TYPE(dft_control_type), POINTER                    :: dft_control
     866             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     867             :       TYPE(neighbor_list_iterator_p_type), &
     868        1400 :          DIMENSION(:), POINTER                           :: nl_iterator
     869             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     870        1400 :          POINTER                                         :: sab_lrc
     871        1400 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     872        1400 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     873        1400 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     874             :       TYPE(se_int_control_type)                          :: se_int_control
     875             :       TYPE(se_taper_type), POINTER                       :: se_taper
     876             :       TYPE(semi_empirical_control_type), POINTER         :: se_control
     877        1400 :       TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
     878             :       TYPE(semi_empirical_type), POINTER                 :: se_kind_a, se_kind_b
     879             : 
     880        1400 :       CALL timeset(routineN, handle)
     881        1400 :       NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper, &
     882        1400 :                efield0, efield1, efield2, atprop)
     883             : 
     884             :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, &
     885             :                       para_env=para_env, sab_lrc=sab_lrc, atomic_kind_set=atomic_kind_set, &
     886        1400 :                       qs_kind_set=qs_kind_set, particle_set=particle_set, atprop=atprop)
     887             : 
     888             :       ! Parameters
     889        1400 :       CALL initialize_se_taper(se_taper, lr_corr=.TRUE.)
     890        1400 :       se_control => dft_control%qs_control%se_control
     891        1400 :       anag = se_control%analytical_gradients
     892        1400 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
     893        1400 :       atener = atprop%energy
     894             : 
     895             :       CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, &
     896             :                                      do_ewald_gks=.FALSE., integral_screening=se_control%integral_screening, &
     897             :                                      shortrange=.FALSE., max_multipole=se_control%max_multipole, &
     898        1400 :                                      pc_coulomb_int=.TRUE.)
     899             : 
     900        1400 :       nspins = dft_control%nspins
     901        1400 :       CPASSERT(ASSOCIATED(matrix_p))
     902        1400 :       CPASSERT(SIZE(ks_matrix) > 0)
     903        1400 :       CPASSERT(ASSOCIATED(diagmat_p))
     904        1400 :       CPASSERT(ASSOCIATED(diagmat_ks))
     905             :       MARK_USED(ks_matrix)
     906             :       MARK_USED(matrix_p)
     907             : 
     908        1400 :       nkind = SIZE(atomic_kind_set)
     909        1400 :       IF (calculate_forces) THEN
     910          20 :          CALL get_qs_env(qs_env=qs_env, force=force)
     911          20 :          delta = se_control%delta
     912          20 :          CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     913             :       END IF
     914             : 
     915             :       ! Allocate arrays for storing partial information on potential, field, field gradient
     916        1400 :       size1 = SIZE(se_nddo_mpole%efield0)
     917        4200 :       ALLOCATE (efield0(size1))
     918       12582 :       efield0 = 0.0_dp
     919        1400 :       size1 = SIZE(se_nddo_mpole%efield1, 1)
     920        1400 :       size2 = SIZE(se_nddo_mpole%efield1, 2)
     921        5600 :       ALLOCATE (efield1(size1, size2))
     922       46128 :       efield1 = 0.0_dp
     923        1400 :       size1 = SIZE(se_nddo_mpole%efield2, 1)
     924        1400 :       size2 = SIZE(se_nddo_mpole%efield2, 2)
     925        5600 :       ALLOCATE (efield2(size1, size2))
     926      113220 :       efield2 = 0.0_dp
     927             : 
     928             :       ! Initialize if virial is requested
     929        1400 :       IF (use_virial) THEN
     930           2 :          ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
     931           2 :          ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
     932           2 :          ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
     933             :       END IF
     934             : 
     935             :       ! Start of the loop for the correction of the pair interactions
     936        1400 :       ecore2 = 0.0_dp
     937        1400 :       enuclear = 0.0_dp
     938        1400 :       itype = get_se_type(dft_control%qs_control%method_id)
     939             : 
     940        9800 :       ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
     941        3848 :       DO ikind = 1, nkind
     942        2448 :          CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
     943        2448 :          se_kind_list(ikind)%se_param => se_kind_a
     944        2448 :          CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
     945        2448 :          se_defined(ikind) = (defined .AND. natorb_a >= 1)
     946        6296 :          se_natorb(ikind) = natorb_a
     947             :       END DO
     948             : 
     949        1400 :       CALL neighbor_list_iterator_create(nl_iterator, sab_lrc)
     950     1432964 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     951     1431564 :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
     952     1431564 :          IF (.NOT. se_defined(ikind)) CYCLE
     953     1431564 :          IF (.NOT. se_defined(jkind)) CYCLE
     954     1431564 :          se_kind_a => se_kind_list(ikind)%se_param
     955     1431564 :          se_kind_b => se_kind_list(jkind)%se_param
     956     1431564 :          natorb_a = se_natorb(ikind)
     957     1431564 :          natorb_b = se_natorb(jkind)
     958     1431564 :          natorb_a2 = natorb_a**2
     959     1431564 :          natorb_b2 = natorb_b**2
     960             : 
     961     1431564 :          IF (inode == 1) THEN
     962             :             CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
     963        9998 :                                    row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
     964        9998 :             CPASSERT(ASSOCIATED(pa_block_a))
     965        9998 :             check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
     966           0 :             CPASSERT(check)
     967             :             CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
     968        9998 :                                    row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
     969        9998 :             CPASSERT(ASSOCIATED(ksa_block_a))
     970      166686 :             p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
     971       19996 :             pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/))
     972        9998 :             IF (nspins == 2) THEN
     973             :                CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
     974           0 :                                       row=iatom, col=iatom, BLOCK=pa_block_b, found=found)
     975           0 :                CPASSERT(ASSOCIATED(pa_block_b))
     976           0 :                check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
     977           0 :                CPASSERT(check)
     978             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
     979           0 :                                       row=iatom, col=iatom, BLOCK=ksa_block_b, found=found)
     980           0 :                CPASSERT(ASSOCIATED(ksa_block_b))
     981           0 :                p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
     982           0 :                pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/))
     983             :             END IF
     984             : 
     985             :          END IF
     986             : 
     987     5726256 :          dr1 = DOT_PRODUCT(rij, rij)
     988     1432964 :          IF (dr1 > rij_threshold) THEN
     989             :             ! Determine the order of the atoms, and in case switch them..
     990     1425973 :             IF (iatom <= jatom) THEN
     991      767729 :                switch = .FALSE.
     992             :             ELSE
     993      658244 :                switch = .TRUE.
     994             :             END IF
     995             : 
     996             :             ! Point-like interaction corrections
     997             :             CALL se_coulomb_ij_interaction(iatom, jatom, task, do_forces=calculate_forces, &
     998             :                                            do_efield=.TRUE., do_stress=use_virial, charges=se_nddo_mpole%charge, &
     999             :                                            dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, &
    1000             :                                            force_ab=force_ab0, efield0=efield0, efield1=efield1, efield2=efield2, &
    1001             :                                            rab2=dr1, rab=rij, ptens11=ptens11, ptens12=ptens12, ptens13=ptens13, &
    1002             :                                            ptens21=ptens21, ptens22=ptens22, ptens23=ptens23, ptens31=ptens31, &
    1003     1425973 :                                            ptens32=ptens32, ptens33=ptens33)
    1004             : 
    1005             :             ! Retrieve blocks for KS and P
    1006             :             CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
    1007     1425973 :                                    row=jatom, col=jatom, BLOCK=pb_block_a, found=found)
    1008     1425973 :             CPASSERT(ASSOCIATED(pb_block_a))
    1009     1425973 :             check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
    1010           0 :             CPASSERT(check)
    1011             :             CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
    1012     1425973 :                                    row=jatom, col=jatom, BLOCK=ksb_block_a, found=found)
    1013     1425973 :             CPASSERT(ASSOCIATED(ksb_block_a))
    1014    23202075 :             p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
    1015     2851946 :             pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/))
    1016             :             ! Handle more than one configuration
    1017     1425973 :             IF (nspins == 2) THEN
    1018             :                CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
    1019           0 :                                       row=jatom, col=jatom, BLOCK=pb_block_b, found=found)
    1020           0 :                CPASSERT(ASSOCIATED(pb_block_b))
    1021           0 :                check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
    1022           0 :                CPASSERT(check)
    1023             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
    1024           0 :                                       row=jatom, col=jatom, BLOCK=ksb_block_b, found=found)
    1025           0 :                CPASSERT(ASSOCIATED(ksb_block_b))
    1026           0 :                check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
    1027           0 :                CPASSERT(check)
    1028           0 :                p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
    1029           0 :                pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/))
    1030             :             END IF
    1031             : 
    1032     2851946 :             SELECT CASE (dft_control%qs_control%method_id)
    1033             :             CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
    1034             :                   do_method_rm1, do_method_mndod)
    1035             :                ! Evaluate nuclear contribution..
    1036             :                CALL corecore_el(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, &
    1037     1425973 :                                 se_int_control=se_int_control, se_taper=se_taper)
    1038     1425973 :                enuclear = enuclear + enuc
    1039             : 
    1040             :                ! Two-centers One-electron terms
    1041     1425973 :                IF (nspins == 1) THEN
    1042     1425973 :                   ecab = 0._dp
    1043             :                   CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
    1044             :                                  pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
    1045     1425973 :                                  se_taper=se_taper, store_int_env=store_int_env)
    1046     1425973 :                   ecore2 = ecore2 + ecab(1) + ecab(2)
    1047           0 :                ELSE IF (nspins == 2) THEN
    1048           0 :                   ecab = 0._dp
    1049             :                   CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
    1050             :                                  pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag, &
    1051           0 :                                  se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
    1052             :                   CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_b, ksb_block_b, &
    1053             :                                  pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
    1054           0 :                                  se_taper=se_taper, store_int_env=store_int_env)
    1055           0 :                   ecore2 = ecore2 + ecab(1) + ecab(2)
    1056             :                END IF
    1057     1425973 :                IF (atener) THEN
    1058           0 :                   atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1) + 0.5_dp*enuc
    1059           0 :                   atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2) + 0.5_dp*enuc
    1060             :                END IF
    1061             :                ! Coulomb Terms
    1062     1425973 :                IF (nspins == 1) THEN
    1063             :                   CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
    1064             :                               ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
    1065     1425973 :                               store_int_env=store_int_env)
    1066           0 :                ELSE IF (nspins == 2) THEN
    1067             :                   CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
    1068             :                               ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
    1069           0 :                               store_int_env=store_int_env)
    1070             : 
    1071             :                   CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, &
    1072             :                               ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
    1073           0 :                               store_int_env=store_int_env)
    1074             :                END IF
    1075             : 
    1076     1425973 :                IF (calculate_forces) THEN
    1077       43876 :                   atom_a = atom_of_kind(iatom)
    1078       43876 :                   atom_b = atom_of_kind(jatom)
    1079             : 
    1080             :                   ! Evaluate nuclear contribution..
    1081             :                   CALL dcorecore_el(se_kind_a, se_kind_b, rij, denuc=force_ab, itype=itype, delta=delta, &
    1082       43876 :                                     anag=anag, se_int_control=se_int_control, se_taper=se_taper)
    1083             : 
    1084             :                   ! Derivatives of the Two-centre One-electron terms
    1085       43876 :                   IF (nspins == 1) THEN
    1086             :                      CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_a, pb_a, itype=itype, anag=anag, &
    1087             :                                      se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
    1088       43876 :                                      delta=delta)
    1089           0 :                   ELSE IF (nspins == 2) THEN
    1090             :                      CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_block_a, pb_block_a, itype=itype, anag=anag, &
    1091           0 :                                      se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
    1092             :                      CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_b, pb_b, itype=itype, anag=anag, &
    1093           0 :                                      se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
    1094             :                   END IF
    1095       43876 :                   IF (use_virial) THEN
    1096        2376 :                      CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
    1097             :                   END IF
    1098      175504 :                   force_ab = force_ab + force_ab0
    1099             : 
    1100             :                   ! Sum up force components
    1101       43876 :                   force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
    1102       43876 :                   force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
    1103             : 
    1104       43876 :                   force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
    1105       43876 :                   force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
    1106             : 
    1107       43876 :                   force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
    1108       43876 :                   force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
    1109             : 
    1110             :                   ! Derivatives of the Coulomb Terms
    1111       43876 :                   force_ab = 0._dp
    1112       43876 :                   IF (nspins == 1) THEN
    1113             :                      CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, &
    1114       43876 :                                   anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
    1115           0 :                   ELSE IF (nspins == 2) THEN
    1116             :                      CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
    1117           0 :                                   anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
    1118             : 
    1119             :                      CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
    1120           0 :                                   anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
    1121             :                   END IF
    1122       43876 :                   IF (switch) THEN
    1123       21260 :                      force_ab(1) = -force_ab(1)
    1124       21260 :                      force_ab(2) = -force_ab(2)
    1125       21260 :                      force_ab(3) = -force_ab(3)
    1126             :                   END IF
    1127       43876 :                   IF (use_virial) THEN
    1128        2376 :                      CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
    1129             :                   END IF
    1130             : 
    1131             :                   ! Sum up force components
    1132       43876 :                   force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
    1133       43876 :                   force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
    1134             : 
    1135       43876 :                   force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
    1136       43876 :                   force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
    1137             : 
    1138       43876 :                   force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
    1139       43876 :                   force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
    1140             :                END IF
    1141             :             CASE DEFAULT
    1142     1425973 :                CPABORT("")
    1143             :             END SELECT
    1144             :          END IF
    1145             :       END DO
    1146        1400 :       CALL neighbor_list_iterator_release(nl_iterator)
    1147             : 
    1148        1400 :       DEALLOCATE (se_kind_list, se_defined, se_natorb)
    1149             : 
    1150             :       ! Sum-up Virial constribution (long-range correction)
    1151        1400 :       IF (use_virial) THEN
    1152           2 :          pv_glob(1, 1) = pv_glob(1, 1) + ptens11
    1153           2 :          pv_glob(1, 2) = pv_glob(1, 2) + (ptens12 + ptens21)*0.5_dp
    1154           2 :          pv_glob(1, 3) = pv_glob(1, 3) + (ptens13 + ptens31)*0.5_dp
    1155           2 :          pv_glob(2, 1) = pv_glob(1, 2)
    1156           2 :          pv_glob(2, 2) = pv_glob(2, 2) + ptens22
    1157           2 :          pv_glob(2, 3) = pv_glob(2, 3) + (ptens23 + ptens32)*0.5_dp
    1158           2 :          pv_glob(3, 1) = pv_glob(1, 3)
    1159           2 :          pv_glob(3, 2) = pv_glob(2, 3)
    1160           2 :          pv_glob(3, 3) = pv_glob(3, 3) + ptens33
    1161             :       END IF
    1162             : 
    1163             :       ! collect information on potential, field, field gradient
    1164       23764 :       CALL para_env%sum(efield0)
    1165       90856 :       CALL para_env%sum(efield1)
    1166      225040 :       CALL para_env%sum(efield2)
    1167       23764 :       se_nddo_mpole%efield0 = se_nddo_mpole%efield0 - efield0
    1168       90856 :       se_nddo_mpole%efield1 = se_nddo_mpole%efield1 - efield1
    1169      225040 :       se_nddo_mpole%efield2 = se_nddo_mpole%efield2 - efield2
    1170             :       ! deallocate working arrays
    1171        1400 :       DEALLOCATE (efield0)
    1172        1400 :       DEALLOCATE (efield1)
    1173        1400 :       DEALLOCATE (efield2)
    1174             : 
    1175             :       ! Corrections for two-centers one-electron terms + nuclear terms
    1176        1400 :       CALL para_env%sum(enuclear)
    1177        1400 :       CALL para_env%sum(ecore2)
    1178        1400 :       energy%hartree = energy%hartree + ecore2
    1179        1400 :       energy%core_overlap = enuclear
    1180        1400 :       CALL finalize_se_taper(se_taper)
    1181        1400 :       CALL timestop(handle)
    1182        2800 :    END SUBROUTINE build_fock_matrix_coul_lrc
    1183             : 
    1184             : ! **************************************************************************************************
    1185             : !> \brief Construction of the residual part (1/R^3) of the Coulomb long-range
    1186             : !>        term of the Fock matrix
    1187             : !>        The 1/R^3 correction works in real-space strictly on the zero-cell,
    1188             : !>        in order to avoid more parameters to be provided in the input..
    1189             : !> \param qs_env ...
    1190             : !> \param ks_matrix ...
    1191             : !> \param matrix_p ...
    1192             : !> \param energy ...
    1193             : !> \param calculate_forces ...
    1194             : !> \author Teodoro Laino [tlaino] - 12.2008
    1195             : ! **************************************************************************************************
    1196           0 :    SUBROUTINE build_fock_matrix_coul_lr_r3(qs_env, ks_matrix, matrix_p, energy, calculate_forces)
    1197             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1198             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix, matrix_p
    1199             :       TYPE(qs_energy_type), POINTER                      :: energy
    1200             :       LOGICAL, INTENT(in)                                :: calculate_forces
    1201             : 
    1202             :       CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_coul_lr_r3'
    1203             : 
    1204             :       INTEGER :: atom_a, atom_b, ewald_type, handle, iatom, ikind, inode, ispin, itype, jatom, &
    1205             :          jkind, natoms, natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nlocal_particles, nspins
    1206           0 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, se_natorb
    1207             :       LOGICAL                                            :: anag, atener, check, defined, found, &
    1208             :                                                             switch
    1209           0 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: se_defined
    1210             :       REAL(KIND=dp)                                      :: dr1, ecore2, r2inv, r3inv, rinv
    1211             :       REAL(KIND=dp), DIMENSION(2)                        :: ecab
    1212             :       REAL(KIND=dp), DIMENSION(2025)                     :: pa_a, pa_b, pb_a, pb_b
    1213             :       REAL(KIND=dp), DIMENSION(3)                        :: dr3inv, force_ab, rij
    1214             :       REAL(KIND=dp), DIMENSION(45, 45)                   :: p_block_tot_a, p_block_tot_b
    1215           0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ksa_block_a, ksa_block_b, ksb_block_a, &
    1216           0 :                                                             ksb_block_b, pa_block_a, pa_block_b, &
    1217           0 :                                                             pb_block_a, pb_block_b
    1218           0 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1219             :       TYPE(atprop_type), POINTER                         :: atprop
    1220             :       TYPE(cell_type), POINTER                           :: cell
    1221             :       TYPE(cp_logger_type), POINTER                      :: logger
    1222           0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: diagmat_ks, diagmat_p
    1223             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1224             :       TYPE(distribution_1d_type), POINTER                :: local_particles
    1225             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
    1226             :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
    1227             :       TYPE(fist_nonbond_env_type), POINTER               :: se_nonbond_env
    1228             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1229             :       TYPE(nddo_mpole_type), POINTER                     :: se_nddo_mpole
    1230             :       TYPE(neighbor_list_iterator_p_type), &
    1231           0 :          DIMENSION(:), POINTER                           :: nl_iterator
    1232             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1233           0 :          POINTER                                         :: sab_orb
    1234           0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1235           0 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
    1236           0 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1237             :       TYPE(section_vals_type), POINTER                   :: se_section
    1238             :       TYPE(semi_empirical_control_type), POINTER         :: se_control
    1239           0 :       TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
    1240             :       TYPE(semi_empirical_type), POINTER                 :: se_kind_a, se_kind_b
    1241             : 
    1242           0 :       CALL timeset(routineN, handle)
    1243             : 
    1244           0 :       CALL get_qs_env(qs_env=qs_env) !sm->dbcsr
    1245             : 
    1246           0 :       NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, &
    1247           0 :                diagmat_p, local_particles, se_control, ewald_env, ewald_pw, &
    1248           0 :                se_nddo_mpole, se_nonbond_env, se_section, sab_orb, logger, atprop)
    1249             : 
    1250           0 :       logger => cp_get_default_logger()
    1251             :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env, &
    1252             :                       atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set, &
    1253             :                       ewald_env=ewald_env, atprop=atprop, &
    1254             :                       local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole, &
    1255           0 :                       se_nonbond_env=se_nonbond_env, sab_orb=sab_orb)
    1256             : 
    1257           0 :       nlocal_particles = SUM(local_particles%n_el(:))
    1258           0 :       natoms = SIZE(particle_set)
    1259           0 :       CALL ewald_env_get(ewald_env, ewald_type=ewald_type)
    1260           0 :       SELECT CASE (ewald_type)
    1261             :       CASE (do_ewald_ewald)
    1262             :          ! Do Nothing
    1263             :       CASE DEFAULT
    1264           0 :          CPABORT("Periodic SE implemented only for standard EWALD sums.")
    1265             :       END SELECT
    1266             : 
    1267             :       ! Parameters
    1268           0 :       se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE")
    1269           0 :       se_control => dft_control%qs_control%se_control
    1270           0 :       anag = se_control%analytical_gradients
    1271           0 :       atener = atprop%energy
    1272             : 
    1273           0 :       nspins = dft_control%nspins
    1274           0 :       CPASSERT(ASSOCIATED(matrix_p))
    1275           0 :       CPASSERT(SIZE(ks_matrix) > 0)
    1276             : 
    1277           0 :       CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
    1278           0 :       CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
    1279             : 
    1280           0 :       nkind = SIZE(atomic_kind_set)
    1281           0 :       DO ispin = 1, nspins
    1282             :          ! Allocate diagonal block matrices
    1283           0 :          ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
    1284           0 :          CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
    1285           0 :          CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
    1286           0 :          CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
    1287           0 :          CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
    1288           0 :          CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
    1289             :       END DO
    1290             : 
    1291             :       ! Possibly compute forces
    1292           0 :       IF (calculate_forces) THEN
    1293           0 :          CALL get_qs_env(qs_env=qs_env, force=force)
    1294           0 :          CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
    1295             :       END IF
    1296             :       itype = get_se_type(dft_control%qs_control%method_id)
    1297             : 
    1298           0 :       ecore2 = 0.0_dp
    1299             :       ! Real space part of the 1/R^3 sum
    1300           0 :       ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
    1301           0 :       DO ikind = 1, nkind
    1302           0 :          CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
    1303           0 :          se_kind_list(ikind)%se_param => se_kind_a
    1304           0 :          CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
    1305           0 :          se_defined(ikind) = (defined .AND. natorb_a >= 1)
    1306           0 :          se_natorb(ikind) = natorb_a
    1307             :       END DO
    1308             : 
    1309           0 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
    1310           0 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    1311           0 :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
    1312           0 :          IF (.NOT. se_defined(ikind)) CYCLE
    1313           0 :          IF (.NOT. se_defined(jkind)) CYCLE
    1314           0 :          se_kind_a => se_kind_list(ikind)%se_param
    1315           0 :          se_kind_b => se_kind_list(jkind)%se_param
    1316           0 :          natorb_a = se_natorb(ikind)
    1317           0 :          natorb_b = se_natorb(jkind)
    1318           0 :          natorb_a2 = natorb_a**2
    1319           0 :          natorb_b2 = natorb_b**2
    1320             : 
    1321           0 :          IF (inode == 1) THEN
    1322             :             CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
    1323           0 :                                    row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
    1324           0 :             CPASSERT(ASSOCIATED(pa_block_a))
    1325           0 :             check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
    1326           0 :             CPASSERT(check)
    1327             :             CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
    1328           0 :                                    row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
    1329           0 :             CPASSERT(ASSOCIATED(ksa_block_a))
    1330           0 :             p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
    1331           0 :             pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/))
    1332           0 :             IF (nspins == 2) THEN
    1333             :                CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
    1334           0 :                                       row=iatom, col=iatom, BLOCK=pa_block_b, found=found)
    1335           0 :                CPASSERT(ASSOCIATED(pa_block_b))
    1336           0 :                check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
    1337           0 :                CPASSERT(check)
    1338             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
    1339           0 :                                       row=iatom, col=iatom, BLOCK=ksa_block_b, found=found)
    1340           0 :                CPASSERT(ASSOCIATED(ksa_block_b))
    1341           0 :                p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
    1342           0 :                pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/))
    1343             :             END IF
    1344             :          END IF
    1345             : 
    1346           0 :          dr1 = DOT_PRODUCT(rij, rij)
    1347           0 :          IF (dr1 > rij_threshold) THEN
    1348             :             ! Determine the order of the atoms, and in case switch them..
    1349           0 :             IF (iatom <= jatom) THEN
    1350           0 :                switch = .FALSE.
    1351             :             ELSE
    1352           0 :                switch = .TRUE.
    1353             :             END IF
    1354             :             ! Retrieve blocks for KS and P
    1355             :             CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
    1356           0 :                                    row=jatom, col=jatom, BLOCK=pb_block_a, found=found)
    1357           0 :             CPASSERT(ASSOCIATED(pb_block_a))
    1358           0 :             check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
    1359           0 :             CPASSERT(check)
    1360             :             CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
    1361           0 :                                    row=jatom, col=jatom, BLOCK=ksb_block_a, found=found)
    1362           0 :             CPASSERT(ASSOCIATED(ksb_block_a))
    1363           0 :             p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
    1364           0 :             pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/))
    1365             :             ! Handle more than one configuration
    1366           0 :             IF (nspins == 2) THEN
    1367             :                CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
    1368           0 :                                       row=jatom, col=jatom, BLOCK=pb_block_b, found=found)
    1369           0 :                CPASSERT(ASSOCIATED(pb_block_b))
    1370           0 :                check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
    1371           0 :                CPASSERT(check)
    1372             :                CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
    1373           0 :                                       row=jatom, col=jatom, BLOCK=ksb_block_b, found=found)
    1374           0 :                CPASSERT(ASSOCIATED(ksb_block_b))
    1375           0 :                check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
    1376           0 :                CPASSERT(check)
    1377           0 :                p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
    1378           0 :                pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/))
    1379             :             END IF
    1380             : 
    1381           0 :             SELECT CASE (dft_control%qs_control%method_id)
    1382             :             CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
    1383             :                   do_method_rm1, do_method_mndod, do_method_pnnl)
    1384             : 
    1385             :                ! Pre-compute some quantities..
    1386           0 :                r2inv = 1.0_dp/dr1
    1387           0 :                rinv = SQRT(r2inv)
    1388           0 :                r3inv = rinv**3
    1389             : 
    1390             :                ! Two-centers One-electron terms
    1391           0 :                IF (nspins == 1) THEN
    1392           0 :                   ecab = 0._dp
    1393             :                   CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_a, pb_a, &
    1394             :                                     ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, &
    1395           0 :                                     e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
    1396           0 :                   ecore2 = ecore2 + ecab(1) + ecab(2)
    1397           0 :                ELSE IF (nspins == 2) THEN
    1398           0 :                   ecab = 0._dp
    1399             :                   CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_block_a, &
    1400             :                                     pb_block_a, ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, &
    1401           0 :                                     e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
    1402             : 
    1403             :                   CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_b, ksb_block_b, pa_b, pb_b, &
    1404             :                                     ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, &
    1405           0 :                                     e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
    1406           0 :                   ecore2 = ecore2 + ecab(1) + ecab(2)
    1407             :                END IF
    1408           0 :                IF (atener) THEN
    1409           0 :                   atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1)
    1410           0 :                   atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2)
    1411             :                END IF
    1412             :                ! Coulomb Terms
    1413           0 :                IF (nspins == 1) THEN
    1414             :                   CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
    1415           0 :                                  ksb_block_a, factor=0.5_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
    1416           0 :                ELSE IF (nspins == 2) THEN
    1417             :                   CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
    1418           0 :                                  ksb_block_a, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
    1419             : 
    1420             :                   CALL fock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, &
    1421           0 :                                  ksb_block_b, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
    1422             :                END IF
    1423             : 
    1424             :                ! Compute forces if requested
    1425           0 :                IF (calculate_forces) THEN
    1426           0 :                   dr3inv = -3.0_dp*rij*r3inv*r2inv
    1427           0 :                   atom_a = atom_of_kind(iatom)
    1428           0 :                   atom_b = atom_of_kind(jatom)
    1429             : 
    1430           0 :                   force_ab = 0.0_dp
    1431             :                   ! Derivatives of the One-centre One-electron terms
    1432           0 :                   IF (nspins == 1) THEN
    1433             :                      CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_a, pb_a, force_ab, &
    1434           0 :                                         se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
    1435           0 :                   ELSE IF (nspins == 2) THEN
    1436             :                      CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_block_a, pb_block_a, force_ab, &
    1437           0 :                                         se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
    1438             : 
    1439             :                      CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_b, pb_b, force_ab, &
    1440           0 :                                         se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
    1441             :                   END IF
    1442             : 
    1443             :                   ! Sum up force components
    1444           0 :                   force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
    1445           0 :                   force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
    1446             : 
    1447           0 :                   force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
    1448           0 :                   force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
    1449             : 
    1450           0 :                   force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
    1451           0 :                   force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
    1452             : 
    1453             :                   ! Derivatives of the Coulomb Terms
    1454           0 :                   force_ab = 0.0_dp
    1455           0 :                   IF (nspins == 1) THEN
    1456             :                      CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, &
    1457           0 :                                      w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
    1458           0 :                   ELSE IF (nspins == 2) THEN
    1459             :                      CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
    1460           0 :                                      w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
    1461             : 
    1462             :                      CALL dfock2C_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
    1463           0 :                                      w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
    1464             :                   END IF
    1465             : 
    1466             :                   ! Sum up force components
    1467           0 :                   force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
    1468           0 :                   force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
    1469             : 
    1470           0 :                   force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
    1471           0 :                   force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
    1472             : 
    1473           0 :                   force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
    1474           0 :                   force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
    1475             :                END IF
    1476             :             CASE DEFAULT
    1477           0 :                CPABORT("")
    1478             :             END SELECT
    1479             :          END IF
    1480             :       END DO
    1481           0 :       CALL neighbor_list_iterator_release(nl_iterator)
    1482             : 
    1483           0 :       DEALLOCATE (se_kind_list, se_defined, se_natorb)
    1484             : 
    1485           0 :       DO ispin = 1, nspins
    1486           0 :          CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
    1487           0 :          CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
    1488             :          CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
    1489           0 :                         1.0_dp, 1.0_dp)
    1490             :       END DO
    1491           0 :       CALL dbcsr_deallocate_matrix_set(diagmat_p)
    1492           0 :       CALL dbcsr_deallocate_matrix_set(diagmat_ks)
    1493             : 
    1494             :       ! Two-centers one-electron terms
    1495           0 :       CALL para_env%sum(ecore2)
    1496           0 :       energy%hartree = energy%hartree + ecore2
    1497             : 
    1498           0 :       CALL timestop(handle)
    1499           0 :    END SUBROUTINE build_fock_matrix_coul_lr_r3
    1500             : 
    1501             : END MODULE se_fock_matrix_coulomb
    1502             : 

Generated by: LCOV version 1.15