LCOV - code coverage report
Current view: top level - src - qs_ks_atom.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 184 189 97.4 %
Date: 2024-03-28 07:31:50 Functions: 3 3 100.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 routines that build the Kohn-Sham matrix  contributions
      10             : !>      coming from local atomic densities
      11             : ! **************************************************************************************************
      12             : MODULE qs_ks_atom
      13             : 
      14             :    USE ao_util,                         ONLY: trace_r_AxB
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16             :                                               get_atomic_kind,&
      17             :                                               get_atomic_kind_set
      18             :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      19             :                                               gto_basis_set_type
      20             :    USE cp_array_utils,                  ONLY: cp_2d_r_p_type
      21             :    USE cp_control_types,                ONLY: dft_control_type
      22             :    USE dbcsr_api,                       ONLY: dbcsr_get_block_p,&
      23             :                                               dbcsr_p_type
      24             :    USE kinds,                           ONLY: dp,&
      25             :                                               int_8
      26             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      27             :                                               kpoint_type
      28             :    USE message_passing,                 ONLY: mp_para_env_type
      29             :    USE qs_environment_types,            ONLY: get_qs_env,&
      30             :                                               qs_environment_type
      31             :    USE qs_force_types,                  ONLY: qs_force_type
      32             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      33             :                                               get_qs_kind_set,&
      34             :                                               qs_kind_type
      35             :    USE qs_neighbor_list_types,          ONLY: get_iterator_task,&
      36             :                                               neighbor_list_iterate,&
      37             :                                               neighbor_list_iterator_create,&
      38             :                                               neighbor_list_iterator_p_type,&
      39             :                                               neighbor_list_iterator_release,&
      40             :                                               neighbor_list_set_p_type,&
      41             :                                               neighbor_list_task_type
      42             :    USE qs_nl_hash_table_types,          ONLY: nl_hash_table_add,&
      43             :                                               nl_hash_table_create,&
      44             :                                               nl_hash_table_get_from_index,&
      45             :                                               nl_hash_table_is_null,&
      46             :                                               nl_hash_table_obj,&
      47             :                                               nl_hash_table_release,&
      48             :                                               nl_hash_table_status
      49             :    USE qs_oce_methods,                  ONLY: prj_gather
      50             :    USE qs_oce_types,                    ONLY: oce_matrix_type
      51             :    USE qs_rho_atom_types,               ONLY: get_rho_atom,&
      52             :                                               rho_atom_coeff,&
      53             :                                               rho_atom_type
      54             :    USE sap_kind_types,                  ONLY: alist_post_align_blk,&
      55             :                                               alist_pre_align_blk,&
      56             :                                               alist_type,&
      57             :                                               get_alist
      58             :    USE util,                            ONLY: get_limit
      59             :    USE virial_methods,                  ONLY: virial_pair_force
      60             :    USE virial_types,                    ONLY: virial_type
      61             : 
      62             : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
      63             : !$                    omp_get_thread_num, &
      64             : !$                    omp_lock_kind, &
      65             : !$                    omp_init_lock, omp_set_lock, &
      66             : !$                    omp_unset_lock, omp_destroy_lock
      67             : 
      68             : #include "./base/base_uses.f90"
      69             : 
      70             :    IMPLICIT NONE
      71             : 
      72             :    PRIVATE
      73             : 
      74             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_atom'
      75             : 
      76             :    PUBLIC :: update_ks_atom
      77             : 
      78             : CONTAINS
      79             : 
      80             : ! **************************************************************************************************
      81             : !> \brief The correction to the KS matrix due to the GAPW local terms to the hartree and
      82             : !>      XC contributions is here added. The correspondig forces contribution are also calculated
      83             : !>      if required. Each sparse-matrix block A-B is corrected by all the atomic contributions
      84             : !>      centered on atoms C for which the triplet A-C-B exists (they are close enough)
      85             : !>      To this end special lists are used
      86             : !> \param qs_env qs environment, for the lists, the contraction coefficients and the
      87             : !>               pre calculated integrals of the potential with the atomic orbitals
      88             : !> \param ksmat KS matrix, sparse matrix
      89             : !> \param pmat density matrix, sparse matrix, needed only for the forces
      90             : !> \param forces switch for the calculation of the forces on atoms
      91             : !> \param tddft switch for TDDFT linear response
      92             : !> \param rho_atom_external ...
      93             : !> \param kind_set_external ...
      94             : !> \param oce_external ...
      95             : !> \param sab_external ...
      96             : !> \param kscale ...
      97             : !> \param kintegral ...
      98             : !> \param kforce ...
      99             : !> \param fscale ...
     100             : !> \par History
     101             : !>      created [MI]
     102             : !>      the loop over the spins is done internally [03-05,MI]
     103             : !>      Rewrite using new OCE matrices [08.02.09, jhu]
     104             : !>      Add OpenMP [Apr 2016, EPCC]
     105             : !>      Allow for external kind_set, rho_atom_set, oce, sab 12.2019 (A. Bussy)
     106             : ! **************************************************************************************************
     107       21182 :    SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, &
     108             :                              kind_set_external, oce_external, sab_external, kscale, &
     109             :                              kintegral, kforce, fscale)
     110             : 
     111             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     112             :       TYPE(dbcsr_p_type), DIMENSION(*), INTENT(INOUT)    :: ksmat, pmat
     113             :       LOGICAL, INTENT(IN)                                :: forces
     114             :       LOGICAL, INTENT(IN), OPTIONAL                      :: tddft
     115             :       TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
     116             :          POINTER                                         :: rho_atom_external
     117             :       TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
     118             :          POINTER                                         :: kind_set_external
     119             :       TYPE(oce_matrix_type), OPTIONAL, POINTER           :: oce_external
     120             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     121             :          OPTIONAL, POINTER                               :: sab_external
     122             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kscale, kintegral, kforce
     123             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN), OPTIONAL  :: fscale
     124             : 
     125             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'update_ks_atom'
     126             : 
     127             :       INTEGER :: bo(2), handle, ia_kind, iac, iat, iatom, ibc, ikind, img, ip, ispin, ja_kind, &
     128             :          jatom, jkind, ka_kind, kac, katom, kbc, kkind, ldCPC, max_gau, max_nsgf, n_cont_a, &
     129             :          n_cont_b, nat, natom, nimages, nkind, nl_table_num_slots, nsoctot, nspins, slot_num
     130             :       INTEGER(KIND=int_8)                                :: nl_table_key
     131       21182 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     132             :       INTEGER, DIMENSION(3)                              :: cell_b
     133       21182 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, list_a, list_b
     134       21182 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     135             :       LOGICAL                                            :: dista, distb, found, is_entry_null, &
     136             :                                                             is_task_valid, my_tddft, paw_atom, &
     137             :                                                             use_virial
     138       21182 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: a_matrix, p_matrix
     139             :       REAL(dp), DIMENSION(3)                             :: rac, rbc
     140             :       REAL(dp), DIMENSION(3, 3)                          :: force_tmp
     141             :       REAL(kind=dp)                                      :: eps_cpc, factor1, factor2
     142       21182 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: C_int_h, C_int_s, coc
     143       21182 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dCPC_h, dCPC_s
     144       21182 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: PC_h, PC_s
     145             :       REAL(KIND=dp), DIMENSION(2)                        :: force_fac
     146             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_virial_thread
     147       21182 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: C_coeff_hh_a, C_coeff_hh_b, &
     148       21182 :                                                             C_coeff_ss_a, C_coeff_ss_b
     149             :       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
     150       21182 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     151       21182 :       TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER        :: mat_h, mat_p
     152             :       TYPE(dft_control_type), POINTER                    :: dft_control
     153       21182 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     154             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     155             :       TYPE(kpoint_type), POINTER                         :: kpoints
     156             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     157             :       TYPE(neighbor_list_iterator_p_type), &
     158       21182 :          DIMENSION(:), POINTER                           :: nl_iterator
     159             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     160       21182 :          POINTER                                         :: my_sab
     161             :       TYPE(neighbor_list_task_type), POINTER             :: next_task, task
     162             :       TYPE(nl_hash_table_obj)                            :: nl_hash_table
     163             :       TYPE(oce_matrix_type), POINTER                     :: my_oce
     164       21182 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     165       21182 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: my_kind_set
     166       21182 :       TYPE(rho_atom_coeff), DIMENSION(:), POINTER        :: int_local_h, int_local_s
     167       21182 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: my_rho_atom
     168             :       TYPE(rho_atom_type), POINTER                       :: rho_at
     169             :       TYPE(virial_type), POINTER                         :: virial
     170             : 
     171       21182 : !$    INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
     172             : !$    INTEGER                                            :: lock_num
     173             : 
     174       21182 :       CALL timeset(routineN, handle)
     175             : 
     176       21182 :       NULLIFY (my_kind_set, atomic_kind_set, force, my_oce, para_env, my_rho_atom, my_sab)
     177       21182 :       NULLIFY (mat_h, mat_p, dft_control)
     178             : 
     179             :       CALL get_qs_env(qs_env=qs_env, &
     180             :                       qs_kind_set=my_kind_set, &
     181             :                       atomic_kind_set=atomic_kind_set, &
     182             :                       force=force, &
     183             :                       oce=my_oce, &
     184             :                       para_env=para_env, &
     185             :                       rho_atom_set=my_rho_atom, &
     186             :                       virial=virial, &
     187             :                       sab_orb=my_sab, &
     188       21182 :                       dft_control=dft_control)
     189             : 
     190       21182 :       nspins = dft_control%nspins
     191       21182 :       nimages = dft_control%nimages
     192             : 
     193       21182 :       factor1 = 1.0_dp
     194       21182 :       factor2 = 1.0_dp
     195             : 
     196             :       !deal with externals
     197       21182 :       my_tddft = .FALSE.
     198       21182 :       IF (PRESENT(tddft)) my_tddft = tddft
     199        5540 :       IF (my_tddft) THEN
     200        2906 :          IF (nspins == 1) factor1 = 2.0_dp
     201        2906 :          CPASSERT(nimages == 1)
     202             :       END IF
     203       21182 :       IF (PRESENT(kscale)) THEN
     204          68 :          factor1 = factor1*kscale
     205          68 :          factor2 = factor2*kscale
     206             :       END IF
     207       21182 :       IF (PRESENT(kintegral)) factor1 = kintegral
     208       21182 :       IF (PRESENT(kforce)) factor2 = kforce
     209       63546 :       force_fac = 1.0_dp
     210       21182 :       IF (PRESENT(fscale)) force_fac(:) = fscale(:)
     211             : 
     212       21182 :       IF (PRESENT(rho_atom_external)) my_rho_atom => rho_atom_external
     213       21182 :       IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
     214       21182 :       IF (PRESENT(oce_external)) my_oce => oce_external
     215       21182 :       IF (PRESENT(sab_external)) my_sab => sab_external
     216             : 
     217             :       ! kpoint images
     218       21182 :       NULLIFY (cell_to_index)
     219       21182 :       IF (nimages > 1) THEN
     220         202 :          CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
     221         202 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     222             :       END IF
     223             : 
     224       21182 :       eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
     225             : 
     226       21182 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
     227       21182 :       CALL get_qs_kind_set(my_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
     228             : 
     229       21182 :       IF (forces) THEN
     230         970 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
     231         970 :          ldCPC = max_gau
     232         970 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     233             :       ELSE
     234             :          use_virial = .FALSE.
     235             :       END IF
     236             : 
     237       21182 :       pv_virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
     238             : 
     239       21182 :       nkind = SIZE(my_kind_set, 1)
     240             :       ! Collect the local potential contributions from all the processors
     241             :       ASSOCIATE (mepos => para_env%mepos, num_pe => para_env%num_pe)
     242       65046 :       DO ikind = 1, nkind
     243       43864 :          NULLIFY (atom_list)
     244       43864 :          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
     245       43864 :          CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
     246       65046 :          IF (paw_atom) THEN
     247             :             ! gather the atomic block integrals in a more compressed format
     248       39878 :             bo = get_limit(nat, num_pe, mepos)
     249       70519 :             DO iat = bo(1), bo(2)
     250       30641 :                iatom = atom_list(iat)
     251      106084 :                DO ispin = 1, nspins
     252             :                   CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
     253       35565 :                                   my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, my_kind_set(ikind))
     254             :                   CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
     255       66206 :                                   my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, my_kind_set(ikind))
     256             :                END DO
     257             :             END DO
     258             :             ! broadcast the CPC arrays to all processors (replicated data)
     259      119634 :             DO ip = 0, num_pe - 1
     260       79756 :                bo = get_limit(nat, num_pe, ip)
     261      180916 :                DO iat = bo(1), bo(2)
     262       61282 :                   iatom = atom_list(iat)
     263      212168 :                   DO ispin = 1, nspins
     264    56455746 :                      CALL para_env%bcast(my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, ip)
     265    56517028 :                      CALL para_env%bcast(my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, ip)
     266             :                   END DO
     267             :                END DO
     268             :             END DO
     269             :          END IF
     270             :       END DO
     271             :       END ASSOCIATE
     272             : 
     273      107410 :       ALLOCATE (basis_set_list(nkind))
     274       65046 :       DO ikind = 1, nkind
     275       43864 :          CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_set_a)
     276       65046 :          IF (ASSOCIATED(basis_set_a)) THEN
     277       43864 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     278             :          ELSE
     279           0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     280             :          END IF
     281             :       END DO
     282             : 
     283             :       ! build the hash table in serial...
     284             :       ! ... creation ...
     285       21182 :       CALL neighbor_list_iterator_create(nl_iterator, my_sab)
     286       21182 :       nl_table_num_slots = natom*natom/2 ! this is probably not optimal, but it seems a good start
     287       21182 :       CALL nl_hash_table_create(nl_hash_table, nmax=nl_table_num_slots)
     288             :       ! ... and population
     289      685422 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0) ! build hash table in serial, so don't pass mepos
     290      664240 :          ALLOCATE (task) ! They must be deallocated before the nl_hash_table is released
     291      664240 :          CALL get_iterator_task(nl_iterator, task) ! build hash table in serial, so don't pass mepos
     292             :          ! tasks with the same key access the same blocks of H & P
     293      664240 :          IF (task%iatom <= task%jatom) THEN
     294      389253 :             nl_table_key = natom*task%iatom + task%jatom
     295             :          ELSE
     296      274987 :             nl_table_key = natom*task%jatom + task%iatom
     297             :          END IF
     298      664240 :          CALL nl_hash_table_add(nl_hash_table, nl_table_key, task)
     299             :       END DO
     300       21182 :       CALL neighbor_list_iterator_release(nl_iterator)
     301             : 
     302             :       ! Get the total number of (possibly empty) entries in the table
     303       21182 :       CALL nl_hash_table_status(nl_hash_table, nmax=nl_table_num_slots)
     304             : 
     305             : !$OMP PARALLEL DEFAULT(NONE)                                         &
     306             : !$OMP           SHARED(nl_table_num_slots, nl_hash_table             &
     307             : !$OMP                 , max_gau, max_nsgf, nspins, forces            &
     308             : !$OMP                 , basis_set_list, nimages, cell_to_index       &
     309             : !$OMP                 , ksmat, pmat, natom, nkind, my_kind_set, my_oce &
     310             : !$OMP                 , my_rho_atom, factor1, factor2, use_virial    &
     311             : !$OMP                 , atom_of_kind, ldCPC, force, locks, force_fac &
     312             : !$OMP                 )                                              &
     313             : !$OMP          PRIVATE( slot_num, is_entry_null, TASK, is_task_valid &
     314             : !$OMP                 , C_int_h, C_int_s, coc, a_matrix, p_matrix    &
     315             : !$OMP                 , mat_h, mat_p, dCPC_h, dCPC_s, PC_h, PC_s     &
     316             : !$OMP                 , int_local_h, int_local_s                     &
     317             : !$OMP                 , ikind, jkind, iatom, jatom, cell_b           &
     318             : !$OMP                 , basis_set_a, basis_set_b, img                &
     319             : !$OMP                 , found, next_task                             &
     320             : !$OMP                 , kkind, paw_atom                              &
     321             : !$OMP                 , iac, alist_ac, kac, n_cont_a, list_a         &
     322             : !$OMP                 , ibc, alist_bc, kbc, n_cont_b, list_b         &
     323             : !$OMP                 , katom, rho_at, nsoctot                       &
     324             : !$OMP                 , C_coeff_hh_a, C_coeff_ss_a, dista, rac       &
     325             : !$OMP                 , C_coeff_hh_b, C_coeff_ss_b, distb, rbc       &
     326             : !$OMP                 , ia_kind, ja_kind, ka_kind, force_tmp         &
     327             : !$OMP                 )                                              &
     328             : !$OMP        REDUCTION(+:pv_virial_thread                            &
     329       21182 : !$OMP                 )
     330             : 
     331             :       ALLOCATE (C_int_h(max_gau*max_nsgf), C_int_s(max_gau*max_nsgf), coc(max_gau*max_gau), &
     332             :                 a_matrix(max_gau, max_gau), p_matrix(max_nsgf, max_nsgf))
     333             : 
     334             :       ALLOCATE (mat_h(nspins), mat_p(nspins))
     335             :       DO ispin = 1, nspins
     336             :          NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
     337             :       END DO
     338             : 
     339             :       IF (forces) THEN
     340             :          ALLOCATE (dCPC_h(max_gau, max_gau), dCPC_s(max_gau, max_gau), &
     341             :                    PC_h(max_nsgf, max_gau, nspins), PC_s(max_nsgf, max_gau, nspins))
     342             : !$OMP SINGLE
     343             : !$       ALLOCATE (locks(natom*nkind))
     344             : !$OMP END SINGLE
     345             : 
     346             : !$OMP DO
     347             : !$       do lock_num = 1, natom*nkind
     348             : !$          call omp_init_lock(locks(lock_num))
     349             : !$       end do
     350             : !$OMP END DO
     351             :       END IF
     352             : 
     353             :       ! Dynamic schedule to take account of the fact that some slots may be empty
     354             :       ! or contain 1 or more tasks
     355             : !$OMP DO SCHEDULE(DYNAMIC,5)
     356             :       DO slot_num = 1, nl_table_num_slots
     357             :          CALL nl_hash_table_is_null(nl_hash_table, slot_num, is_entry_null)
     358             : 
     359             :          IF (.NOT. is_entry_null) THEN
     360             :             CALL nl_hash_table_get_from_index(nl_hash_table, slot_num, task)
     361             : 
     362             :             is_task_valid = ASSOCIATED(task)
     363             :             DO WHILE (is_task_valid)
     364             : 
     365             :                ikind = task%ikind
     366             :                jkind = task%jkind
     367             :                iatom = task%iatom
     368             :                jatom = task%jatom
     369             :                cell_b(:) = task%cell(:)
     370             : 
     371             :                basis_set_a => basis_set_list(ikind)%gto_basis_set
     372             :                IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     373             :                basis_set_b => basis_set_list(jkind)%gto_basis_set
     374             :                IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     375             : 
     376             :                IF (nimages > 1) THEN
     377             :                   img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
     378             :                   CPASSERT(img > 0)
     379             :                ELSE
     380             :                   img = 1
     381             :                END IF
     382             : 
     383             :                DO ispin = 1, nspins
     384             :                   NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
     385             : 
     386             :                   found = .FALSE.
     387             :                   IF (iatom <= jatom) THEN
     388             :                      CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
     389             :                                             row=iatom, col=jatom, &
     390             :                                             BLOCK=mat_h(ispin)%array, found=found)
     391             :                   ELSE
     392             :                      CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
     393             :                                             row=jatom, col=iatom, &
     394             :                                             BLOCK=mat_h(ispin)%array, found=found)
     395             :                   END IF
     396             :                   CPASSERT(found)
     397             : 
     398             :                   IF (forces) THEN
     399             :                      found = .FALSE.
     400             :                      IF (iatom <= jatom) THEN
     401             :                         CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
     402             :                                                row=iatom, col=jatom, &
     403             :                                                BLOCK=mat_p(ispin)%array, found=found)
     404             :                      ELSE
     405             :                         CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
     406             :                                                row=jatom, col=iatom, &
     407             :                                                BLOCK=mat_p(ispin)%array, found=found)
     408             :                      END IF
     409             :                      CPASSERT(found)
     410             :                   END IF
     411             :                END DO
     412             : 
     413             :                DO kkind = 1, nkind
     414             :                   CALL get_qs_kind(my_kind_set(kkind), paw_atom=paw_atom)
     415             : 
     416             :                   IF (.NOT. paw_atom) CYCLE
     417             : 
     418             :                   iac = ikind + nkind*(kkind - 1)
     419             :                   ibc = jkind + nkind*(kkind - 1)
     420             :                   IF (.NOT. ASSOCIATED(my_oce%intac(iac)%alist)) CYCLE
     421             :                   IF (.NOT. ASSOCIATED(my_oce%intac(ibc)%alist)) CYCLE
     422             : 
     423             :                   CALL get_alist(my_oce%intac(iac), alist_ac, iatom)
     424             :                   CALL get_alist(my_oce%intac(ibc), alist_bc, jatom)
     425             :                   IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
     426             :                   IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
     427             : 
     428             :                   DO kac = 1, alist_ac%nclist
     429             :                      DO kbc = 1, alist_bc%nclist
     430             :                         IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
     431             : 
     432             :                         IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
     433             :                            n_cont_a = alist_ac%clist(kac)%nsgf_cnt
     434             :                            n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
     435             :                            IF (n_cont_a .EQ. 0 .OR. n_cont_b .EQ. 0) CYCLE
     436             : 
     437             :                            list_a => alist_ac%clist(kac)%sgf_list
     438             :                            list_b => alist_bc%clist(kbc)%sgf_list
     439             : 
     440             :                            katom = alist_ac%clist(kac)%catom
     441             : 
     442             :                            IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
     443             :                               C_coeff_hh_a => alist_ac%clist(kac)%achint
     444             :                               C_coeff_ss_a => alist_ac%clist(kac)%acint
     445             :                               dista = .FALSE.
     446             :                            ELSE
     447             :                               C_coeff_hh_a => alist_ac%clist(kac)%acint
     448             :                               C_coeff_ss_a => alist_ac%clist(kac)%acint
     449             :                               dista = .TRUE.
     450             :                            END IF
     451             : 
     452             :                            IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
     453             :                               C_coeff_hh_b => alist_bc%clist(kbc)%achint
     454             :                               C_coeff_ss_b => alist_bc%clist(kbc)%acint
     455             :                               distb = .FALSE.
     456             :                            ELSE
     457             :                               C_coeff_hh_b => alist_bc%clist(kbc)%acint
     458             :                               C_coeff_ss_b => alist_bc%clist(kbc)%acint
     459             :                               distb = .TRUE.
     460             :                            END IF
     461             : 
     462             :                            rho_at => my_rho_atom(katom)
     463             :                            nsoctot = SIZE(C_coeff_ss_a, 2)
     464             : 
     465             :                            CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
     466             :                            CALL add_vhxca_to_ks(mat_h, C_coeff_hh_a, C_coeff_hh_b, C_coeff_ss_a, C_coeff_ss_b, &
     467             :                                                 int_local_h, int_local_s, nspins, iatom, jatom, nsoctot, factor1, &
     468             :                                                 list_a, n_cont_a, list_b, n_cont_b, C_int_h, C_int_s, a_matrix, dista, distb, coc)
     469             : 
     470             :                            IF (forces) THEN
     471             :                               IF (use_virial) THEN
     472             :                                  rac = alist_ac%clist(kac)%rac
     473             :                                  rbc = alist_bc%clist(kbc)%rac
     474             :                               END IF
     475             :                               ia_kind = atom_of_kind(iatom)
     476             :                               ja_kind = atom_of_kind(jatom)
     477             :                               ka_kind = atom_of_kind(katom)
     478             :                               rho_at => my_rho_atom(katom)
     479             :                               force_tmp(1:3, 1:3) = 0.0_dp
     480             : 
     481             :                               CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
     482             :                               IF (iatom <= jatom) THEN
     483             :                                  CALL add_vhxca_forces(mat_p, C_coeff_hh_a, C_coeff_hh_b, C_coeff_ss_a, C_coeff_ss_b, &
     484             :                                                        int_local_h, int_local_s, force_tmp, nspins, iatom, jatom, nsoctot, &
     485             :                                                        list_a, n_cont_a, list_b, n_cont_b, dCPC_h, dCPC_s, ldCPC, &
     486             :                                                        PC_h, PC_s, p_matrix, force_fac)
     487             :                                  force_tmp = factor2*force_tmp
     488             : !$                               CALL omp_set_lock(locks((ka_kind - 1)*nkind + kkind))
     489             :                                  force(kkind)%vhxc_atom(1:3, ka_kind) = force(kkind)%vhxc_atom(1:3, ka_kind) + force_tmp(1:3, 3)
     490             : !$                               CALL omp_unset_lock(locks((ka_kind - 1)*nkind + kkind))
     491             : !$                               CALL omp_set_lock(locks((ia_kind - 1)*nkind + ikind))
     492             :                                  force(ikind)%vhxc_atom(1:3, ia_kind) = force(ikind)%vhxc_atom(1:3, ia_kind) + force_tmp(1:3, 1)
     493             : !$                               CALL omp_unset_lock(locks((ia_kind - 1)*nkind + ikind))
     494             : !$                               CALL omp_set_lock(locks((ja_kind - 1)*nkind + jkind))
     495             :                                  force(jkind)%vhxc_atom(1:3, ja_kind) = force(jkind)%vhxc_atom(1:3, ja_kind) + force_tmp(1:3, 2)
     496             : !$                               CALL omp_unset_lock(locks((ja_kind - 1)*nkind + jkind))
     497             :                                  IF (use_virial) THEN
     498             :                                     CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 1), rac)
     499             :                                     CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 2), rbc)
     500             :                                  END IF
     501             :                               ELSE
     502             :                                  CALL add_vhxca_forces(mat_p, C_coeff_hh_b, C_coeff_hh_a, C_coeff_ss_b, C_coeff_ss_a, &
     503             :                                                        int_local_h, int_local_s, force_tmp, nspins, jatom, iatom, nsoctot, &
     504             :                                                        list_b, n_cont_b, list_a, n_cont_a, dCPC_h, dCPC_s, ldCPC, &
     505             :                                                        PC_h, PC_s, p_matrix, force_fac)
     506             :                                  force_tmp = factor2*force_tmp
     507             : !$                               CALL omp_set_lock(locks((ka_kind - 1)*nkind + kkind))
     508             :                                  force(kkind)%vhxc_atom(1:3, ka_kind) = force(kkind)%vhxc_atom(1:3, ka_kind) + force_tmp(1:3, 3)
     509             : !$                               CALL omp_unset_lock(locks((ka_kind - 1)*nkind + kkind))
     510             : !$                               CALL omp_set_lock(locks((ia_kind - 1)*nkind + ikind))
     511             :                                  force(ikind)%vhxc_atom(1:3, ia_kind) = force(ikind)%vhxc_atom(1:3, ia_kind) + force_tmp(1:3, 2)
     512             : !$                               CALL omp_unset_lock(locks((ia_kind - 1)*nkind + ikind))
     513             : !$                               CALL omp_set_lock(locks((ja_kind - 1)*nkind + jkind))
     514             :                                  force(jkind)%vhxc_atom(1:3, ja_kind) = force(jkind)%vhxc_atom(1:3, ja_kind) + force_tmp(1:3, 1)
     515             : !$                               CALL omp_unset_lock(locks((ja_kind - 1)*nkind + jkind))
     516             :                                  IF (use_virial) THEN
     517             :                                     CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 2), rac)
     518             :                                     CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 1), rbc)
     519             :                                  END IF
     520             :                               END IF
     521             : 
     522             :                            END IF
     523             :                            EXIT ! search loop over jatom-katom list
     524             :                         END IF
     525             :                      END DO ! kbc
     526             :                   END DO ! kac
     527             : 
     528             :                END DO ! kkind
     529             : 
     530             :                next_task => task%next
     531             :                ! We are done with this task, we can deallocate it
     532             :                DEALLOCATE (task)
     533             :                is_task_valid = ASSOCIATED(next_task)
     534             :                IF (is_task_valid) task => next_task
     535             : 
     536             :             END DO
     537             : 
     538             :          ELSE
     539             :             ! NO KEY/VALUE
     540             :          END IF
     541             :       END DO
     542             : !$OMP END DO
     543             : 
     544             :       DO ispin = 1, nspins
     545             :          NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
     546             :       END DO
     547             :       DEALLOCATE (mat_h, mat_p, C_int_h, C_int_s, a_matrix, p_matrix, coc)
     548             : 
     549             :       IF (forces) THEN
     550             :          DEALLOCATE (dCPC_h, dCPC_s, PC_h, PC_s)
     551             : 
     552             :          ! Implicit barrier at end of OMP DO, so locks can be freed
     553             : !$OMP DO
     554             : !$       DO lock_num = 1, natom*nkind
     555             : !$          call omp_destroy_lock(locks(lock_num))
     556             : !$       END DO
     557             : !$OMP END DO
     558             : 
     559             : !$OMP SINGLE
     560             : !$       DEALLOCATE (locks)
     561             : !$OMP END SINGLE NOWAIT
     562             :       END IF
     563             : 
     564             : !$OMP END PARALLEL
     565             : 
     566       21182 :       IF (use_virial) THEN
     567         364 :          virial%pv_gapw(1:3, 1:3) = virial%pv_gapw(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
     568         364 :          virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
     569             :       END IF
     570             : 
     571       21182 :       CALL nl_hash_table_release(nl_hash_table)
     572             : 
     573       21182 :       DEALLOCATE (basis_set_list)
     574             : 
     575       21182 :       CALL timestop(handle)
     576             : 
     577       84728 :    END SUBROUTINE update_ks_atom
     578             : 
     579             : ! **************************************************************************************************
     580             : !> \brief ...
     581             : !> \param mat_h ...
     582             : !> \param C_hh_a ...
     583             : !> \param C_hh_b ...
     584             : !> \param C_ss_a ...
     585             : !> \param C_ss_b ...
     586             : !> \param int_local_h ...
     587             : !> \param int_local_s ...
     588             : !> \param nspins ...
     589             : !> \param ia ...
     590             : !> \param ja ...
     591             : !> \param nsp ...
     592             : !> \param factor ...
     593             : !> \param lista ...
     594             : !> \param nconta ...
     595             : !> \param listb ...
     596             : !> \param ncontb ...
     597             : !> \param C_int_h ...
     598             : !> \param C_int_s ...
     599             : !> \param a_matrix ...
     600             : !> \param dista ...
     601             : !> \param distb ...
     602             : !> \param coc ...
     603             : ! **************************************************************************************************
     604     4292602 :    SUBROUTINE add_vhxca_to_ks(mat_h, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
     605             :                               int_local_h, int_local_s, &
     606     4292602 :                               nspins, ia, ja, nsp, factor, lista, nconta, listb, ncontb, &
     607     8585204 :                               C_int_h, C_int_s, a_matrix, dista, distb, coc)
     608             :       TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER        :: mat_h
     609             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: C_hh_a, C_hh_b, C_ss_a, C_ss_b
     610             :       TYPE(rho_atom_coeff), DIMENSION(:), POINTER        :: int_local_h, int_local_s
     611             :       INTEGER, INTENT(IN)                                :: nspins, ia, ja, nsp
     612             :       REAL(KIND=dp), INTENT(IN)                          :: factor
     613             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lista
     614             :       INTEGER, INTENT(IN)                                :: nconta
     615             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: listb
     616             :       INTEGER, INTENT(IN)                                :: ncontb
     617             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: C_int_h, C_int_s
     618             :       REAL(dp), DIMENSION(:, :)                          :: a_matrix
     619             :       LOGICAL, INTENT(IN)                                :: dista, distb
     620             :       REAL(dp), DIMENSION(:)                             :: coc
     621             : 
     622             :       INTEGER                                            :: i, ispin, j, k
     623     4292602 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: h_block, int_hard, int_soft
     624             : 
     625     4292602 :       NULLIFY (int_hard, int_soft)
     626             : 
     627     8766757 :       DO ispin = 1, nspins
     628     4474155 :          h_block => mat_h(ispin)%array
     629             :          !
     630     4474155 :          int_hard => int_local_h(ispin)%r_coef
     631     4474155 :          int_soft => int_local_s(ispin)%r_coef
     632             :          !
     633     8766757 :          IF (ia <= ja) THEN
     634     2587004 :             IF (dista .AND. distb) THEN
     635     2453249 :                k = 0
     636    53565157 :                DO i = 1, nsp
     637  1325142553 :                   DO j = 1, nsp
     638  1271577396 :                      k = k + 1
     639  1322689304 :                      coc(k) = int_hard(j, i) - int_soft(j, i)
     640             :                   END DO
     641             :                END DO
     642             :                CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, coc, nsp, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     643     2453249 :                           0.0_dp, C_int_h, nsp)
     644             :                CALL DGEMM('N', 'N', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
     645     2453249 :                           C_int_h, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
     646      133755 :             ELSEIF (dista) THEN
     647             :                CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
     648       44683 :                           C_hh_b(:, :, 1), SIZE(C_hh_b, 1), 0.0_dp, coc, nsp)
     649             :                CALL DGEMM('N', 'T', nsp, ncontb, nsp, -1.0_dp, int_soft, SIZE(int_soft, 1), &
     650       44683 :                           C_ss_b(:, :, 1), SIZE(C_ss_b, 1), 1.0_dp, coc, nsp)
     651             :                CALL DGEMM('N', 'N', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
     652       44683 :                           coc, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
     653       89072 :             ELSEIF (distb) THEN
     654             :                CALL DGEMM('N', 'N', nconta, nsp, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
     655       53507 :                           int_hard, SIZE(int_hard, 1), 0.0_dp, coc, nconta)
     656             :                CALL DGEMM('N', 'N', nconta, nsp, nsp, -factor, C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
     657       53507 :                           int_soft, SIZE(int_soft, 1), 1.0_dp, coc, nconta)
     658             :                CALL DGEMM('N', 'T', nconta, ncontb, nsp, 1.0_dp, coc, nconta, &
     659       53507 :                           C_hh_b(:, :, 1), SIZE(C_hh_b, 1), 0.0_dp, a_matrix, SIZE(a_matrix, 1))
     660             :             ELSE
     661             :                CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
     662             :                           C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     663       35565 :                           0.0_dp, C_int_h, nsp)
     664             :                CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, int_soft, SIZE(int_soft, 1), &
     665             :                           C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
     666       35565 :                           0.0_dp, C_int_s, nsp)
     667             :                CALL DGEMM('N', 'N', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
     668             :                           C_int_h, nsp, &
     669       35565 :                           0.0_dp, a_matrix, SIZE(a_matrix, 1))
     670             :                CALL DGEMM('N', 'N', nconta, ncontb, nsp, -factor, C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
     671             :                           C_int_s, nsp, &
     672       35565 :                           1.0_dp, a_matrix, SIZE(a_matrix, 1))
     673             :             END IF
     674             :             !
     675             :             CALL alist_post_align_blk(a_matrix, SIZE(a_matrix, 1), h_block, SIZE(h_block, 1), &
     676     2587004 :                                       lista, nconta, listb, ncontb)
     677             :          ELSE
     678     1887151 :             IF (dista .AND. distb) THEN
     679     1772170 :                k = 0
     680    38218659 :                DO i = 1, nsp
     681   904333266 :                   DO j = 1, nsp
     682   866114607 :                      k = k + 1
     683   902561096 :                      coc(k) = int_hard(j, i) - int_soft(j, i)
     684             :                   END DO
     685             :                END DO
     686     1772170 :                CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, coc, nsp, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), 0.0_dp, C_int_h, nsp)
     687             :                CALL DGEMM('N', 'N', ncontb, nconta, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     688     1772170 :                           C_int_h, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
     689      114981 :             ELSEIF (dista) THEN
     690             :                CALL DGEMM('N', 'N', ncontb, nsp, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     691       63076 :                           int_hard, SIZE(int_hard, 1), 0.0_dp, coc, ncontb)
     692             :                CALL DGEMM('N', 'N', ncontb, nsp, nsp, -factor, C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
     693       63076 :                           int_soft, SIZE(int_soft, 1), 1.0_dp, coc, ncontb)
     694             :                CALL DGEMM('N', 'T', ncontb, nconta, nsp, 1.0_dp, coc, ncontb, &
     695       63076 :                           C_hh_a(:, :, 1), SIZE(C_hh_a, 1), 0.0_dp, a_matrix, SIZE(a_matrix, 1))
     696       51905 :             ELSEIF (distb) THEN
     697             :                CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
     698       51905 :                           C_hh_a(:, :, 1), SIZE(C_hh_a, 1), 0.0_dp, coc, nsp)
     699             :                CALL DGEMM('N', 'T', nsp, nconta, nsp, -1.0_dp, int_soft, SIZE(int_soft, 1), &
     700       51905 :                           C_ss_a(:, :, 1), SIZE(C_ss_a, 1), 1.0_dp, coc, nsp)
     701             :                CALL DGEMM('N', 'N', ncontb, nconta, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     702       51905 :                           coc, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
     703             :             ELSE
     704             :                CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
     705             :                           C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
     706           0 :                           0.0_dp, C_int_h, nsp)
     707             :                CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, int_soft, SIZE(int_soft, 1), &
     708             :                           C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
     709           0 :                           0.0_dp, C_int_s, nsp)
     710             :                CALL DGEMM('N', 'N', ncontb, nconta, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     711             :                           C_int_h, nsp, &
     712           0 :                           0.0_dp, a_matrix, SIZE(a_matrix, 1))
     713             :                CALL DGEMM('N', 'N', ncontb, nconta, nsp, -factor, C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
     714             :                           C_int_s, nsp, &
     715           0 :                           1.0_dp, a_matrix, SIZE(a_matrix, 1))
     716             :             END IF
     717             :             !
     718             :             CALL alist_post_align_blk(a_matrix, SIZE(a_matrix, 1), h_block, SIZE(h_block, 1), &
     719     1887151 :                                       listb, ncontb, lista, nconta)
     720             :          END IF
     721             :       END DO
     722             : 
     723     4292602 :    END SUBROUTINE add_vhxca_to_ks
     724             : 
     725             : ! **************************************************************************************************
     726             : !> \brief ...
     727             : !> \param mat_p ...
     728             : !> \param C_hh_a ...
     729             : !> \param C_hh_b ...
     730             : !> \param C_ss_a ...
     731             : !> \param C_ss_b ...
     732             : !> \param int_local_h ...
     733             : !> \param int_local_s ...
     734             : !> \param force ...
     735             : !> \param nspins ...
     736             : !> \param ia ...
     737             : !> \param ja ...
     738             : !> \param nsp ...
     739             : !> \param lista ...
     740             : !> \param nconta ...
     741             : !> \param listb ...
     742             : !> \param ncontb ...
     743             : !> \param dCPC_h ...
     744             : !> \param dCPC_s ...
     745             : !> \param ldCPC ...
     746             : !> \param PC_h ...
     747             : !> \param PC_s ...
     748             : !> \param p_matrix ...
     749             : !> \param force_scaling ...
     750             : ! **************************************************************************************************
     751     2966445 :    SUBROUTINE add_vhxca_forces(mat_p, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
     752             :                                int_local_h, int_local_s, &
     753      329605 :                                force, nspins, ia, ja, nsp, lista, nconta, listb, ncontb, &
     754      659210 :                                dCPC_h, dCPC_s, ldCPC, PC_h, PC_s, p_matrix, force_scaling)
     755             :       TYPE(cp_2d_r_p_type), DIMENSION(:), INTENT(IN), &
     756             :          POINTER                                         :: mat_p
     757             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: C_hh_a, C_hh_b, C_ss_a, C_ss_b
     758             :       TYPE(rho_atom_coeff), DIMENSION(:), POINTER        :: int_local_h, int_local_s
     759             :       REAL(dp), DIMENSION(3, 3), INTENT(INOUT)           :: force
     760             :       INTEGER, INTENT(IN)                                :: nspins, ia, ja, nsp
     761             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: lista
     762             :       INTEGER, INTENT(IN)                                :: nconta
     763             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: listb
     764             :       INTEGER, INTENT(IN)                                :: ncontb
     765             :       REAL(KIND=dp), DIMENSION(:, :)                     :: dCPC_h, dCPC_s
     766             :       INTEGER, INTENT(IN)                                :: ldCPC
     767             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: PC_h, PC_s
     768             :       REAL(KIND=dp), DIMENSION(:, :)                     :: p_matrix
     769             :       REAL(KIND=dp), DIMENSION(2), INTENT(IN)            :: force_scaling
     770             : 
     771             :       INTEGER                                            :: dir, ispin
     772      329605 :       REAL(dp), DIMENSION(:, :), POINTER                 :: int_hard, int_soft
     773             :       REAL(KIND=dp)                                      :: ieqj, trace
     774             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block
     775             : 
     776             : !   if(dista.and.distb) we could also here use ChPCh = CsPCs
     777             : !   *** factor 2 because only half of the pairs with ia =/ ja are considered
     778             : 
     779      329605 :       ieqj = 2.0_dp
     780       91042 :       IF (ia == ja) ieqj = 1.0_dp
     781             : 
     782      329605 :       NULLIFY (int_hard, int_soft)
     783             : 
     784      661716 :       DO ispin = 1, nspins
     785      332111 :          p_block => mat_p(ispin)%array
     786             : 
     787             :          CALL alist_pre_align_blk(p_block, SIZE(p_block, 1), p_matrix, SIZE(p_matrix, 1), &
     788      332111 :                                   lista, nconta, listb, ncontb)
     789             : 
     790      332111 :          int_hard => int_local_h(ispin)%r_coef
     791      332111 :          int_soft => int_local_s(ispin)%r_coef
     792             : 
     793             :          CALL DGEMM('N', 'N', nconta, nsp, ncontb, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
     794             :                     C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
     795      332111 :                     0.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1))
     796             :          CALL DGEMM('N', 'N', nconta, nsp, ncontb, 1.0_dp, p_matrix(:, :), SIZE(p_matrix, 1), &
     797             :                     C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
     798      332111 :                     0.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1))
     799             : 
     800     1328444 :          DO dir = 2, 4
     801             :             CALL DGEMM('T', 'N', nsp, nsp, nconta, 1.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1), &
     802             :                        C_hh_a(:, :, dir), SIZE(C_hh_a, 1), &
     803      996333 :                        0.0_dp, dCPC_h, SIZE(dCPC_h, 1))
     804      996333 :             trace = trace_r_AxB(dCPC_h, ldCPC, int_hard, nsp, nsp, nsp)
     805      996333 :             force(dir - 1, 3) = force(dir - 1, 3) + ieqj*trace*force_scaling(ispin)
     806      996333 :             force(dir - 1, 1) = force(dir - 1, 1) - ieqj*trace*force_scaling(ispin)
     807             : 
     808             :             CALL DGEMM('T', 'N', nsp, nsp, nconta, 1.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1), &
     809             :                        C_ss_a(:, :, dir), SIZE(C_ss_a, 1), &
     810      996333 :                        0.0_dp, dCPC_s, SIZE(dCPC_s, 1))
     811      996333 :             trace = trace_r_AxB(dCPC_s, ldCPC, int_soft, nsp, nsp, nsp)
     812      996333 :             force(dir - 1, 3) = force(dir - 1, 3) - ieqj*trace*force_scaling(ispin)
     813     1328444 :             force(dir - 1, 1) = force(dir - 1, 1) + ieqj*trace*force_scaling(ispin)
     814             :          END DO
     815             : 
     816             :          ! j-k contributions
     817             :          CALL DGEMM('T', 'N', ncontb, nsp, nconta, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
     818             :                     C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
     819      332111 :                     0.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1))
     820             :          CALL DGEMM('T', 'N', ncontb, nsp, nconta, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
     821             :                     C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
     822      332111 :                     0.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1))
     823             : 
     824     1658049 :          DO dir = 2, 4
     825             :             CALL DGEMM('T', 'N', nsp, nsp, ncontb, 1.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1), &
     826             :                        C_hh_b(:, :, dir), SIZE(C_hh_b, 1), &
     827      996333 :                        0.0_dp, dCPC_h, SIZE(dCPC_h, 1))
     828      996333 :             trace = trace_r_AxB(dCPC_h, ldCPC, int_hard, nsp, nsp, nsp)
     829      996333 :             force(dir - 1, 3) = force(dir - 1, 3) + ieqj*trace*force_scaling(ispin)
     830      996333 :             force(dir - 1, 2) = force(dir - 1, 2) - ieqj*trace*force_scaling(ispin)
     831             : 
     832             :             CALL DGEMM('T', 'N', nsp, nsp, ncontb, 1.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1), &
     833             :                        C_ss_b(:, :, dir), SIZE(C_ss_b, 1), &
     834      996333 :                        0.0_dp, dCPC_s, SIZE(dCPC_s, 1))
     835      996333 :             trace = trace_r_AxB(dCPC_s, ldCPC, int_soft, nsp, nsp, nsp)
     836      996333 :             force(dir - 1, 3) = force(dir - 1, 3) - ieqj*trace*force_scaling(ispin)
     837     1328444 :             force(dir - 1, 2) = force(dir - 1, 2) + ieqj*trace*force_scaling(ispin)
     838             :          END DO
     839             : 
     840             :       END DO !ispin
     841             : 
     842      329605 :    END SUBROUTINE add_vhxca_forces
     843             : 
     844             : END MODULE qs_ks_atom

Generated by: LCOV version 1.15