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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : ! **************************************************************************************************
       8              : !> \brief Calculation of the local potential contribution of the nonadditive kinetic energy
       9              : !>         <a|V(local)|b> = <a|Sum e^a*rc**2|b>
      10              : !> \par History
      11              : !>      - adapted from core_ppl
      12              : ! **************************************************************************************************
      13              : MODULE kg_tnadd_mat
      14              :    USE ai_overlap_ppl,                  ONLY: ppl_integral
      15              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16              :                                               get_atomic_kind_set
      17              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      18              :                                               gto_basis_set_type
      19              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      20              :                                               dbcsr_create,&
      21              :                                               dbcsr_distribution_type,&
      22              :                                               dbcsr_finalize,&
      23              :                                               dbcsr_get_block_p,&
      24              :                                               dbcsr_p_type,&
      25              :                                               dbcsr_set,&
      26              :                                               dbcsr_type_symmetric
      27              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      28              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      29              :                                               dbcsr_deallocate_matrix_set
      30              :    USE external_potential_types,        ONLY: get_potential,&
      31              :                                               local_potential_type
      32              :    USE kg_environment_types,            ONLY: kg_environment_type
      33              :    USE kinds,                           ONLY: dp
      34              :    USE orbital_pointers,                ONLY: init_orbital_pointers,&
      35              :                                               ncoset
      36              :    USE particle_methods,                ONLY: get_particle_set
      37              :    USE particle_types,                  ONLY: particle_type
      38              :    USE qs_force_types,                  ONLY: qs_force_type
      39              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      40              :                                               get_qs_kind_set,&
      41              :                                               qs_kind_type
      42              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      43              :                                               neighbor_list_iterate,&
      44              :                                               neighbor_list_iterator_create,&
      45              :                                               neighbor_list_iterator_p_type,&
      46              :                                               neighbor_list_iterator_release,&
      47              :                                               neighbor_list_set_p_type,&
      48              :                                               nl_set_sub_iterator,&
      49              :                                               nl_sub_iterate
      50              :    USE virial_methods,                  ONLY: virial_pair_force
      51              :    USE virial_types,                    ONLY: virial_type
      52              : 
      53              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      54              : 
      55              : #include "./base/base_uses.f90"
      56              : 
      57              :    IMPLICIT NONE
      58              : 
      59              :    PRIVATE
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_tnadd_mat'
      62              : 
      63              :    PUBLIC :: build_tnadd_mat
      64              : 
      65              : CONTAINS
      66              : 
      67              : !==========================================================================================================
      68              : 
      69              : ! **************************************************************************************************
      70              : !> \brief ...
      71              : !> \param kg_env ...
      72              : !> \param matrix_p ...
      73              : !> \param force ...
      74              : !> \param virial ...
      75              : !> \param calculate_forces ...
      76              : !> \param use_virial ...
      77              : !> \param qs_kind_set ...
      78              : !> \param atomic_kind_set ...
      79              : !> \param particle_set ...
      80              : !> \param sab_orb ...
      81              : !> \param dbcsr_dist ...
      82              : ! **************************************************************************************************
      83           42 :    SUBROUTINE build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
      84              :                               qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
      85              : 
      86              :       TYPE(kg_environment_type), POINTER                 :: kg_env
      87              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p
      88              :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      89              :       TYPE(virial_type), POINTER                         :: virial
      90              :       LOGICAL, INTENT(IN)                                :: calculate_forces
      91              :       LOGICAL                                            :: use_virial
      92              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      93              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      94              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      95              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      96              :          POINTER                                         :: sab_orb
      97              :       TYPE(dbcsr_distribution_type), POINTER             :: dbcsr_dist
      98              : 
      99              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_tnadd_mat'
     100              :       INTEGER, PARAMETER                                 :: nexp_max = 10
     101              : 
     102              :       INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, icol, ikind, img, imol, inode, irow, &
     103              :          iset, jatom, jkind, jmol, jset, katom, kkind, kmol, last_iatom, last_jatom, ldai, ldsab, &
     104              :          maxco, maxder, maxl, maxlgto, maxnset, maxpol, maxsgf, mepos, natom, ncoa, ncob, nder, &
     105              :          ngau, nkind, npol, nseta, nsetb, nthread, sgfa, sgfb
     106           42 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     107           42 :       INTEGER, DIMENSION(:), POINTER                     :: atom_to_molecule, col_blk_sizes, la_max, &
     108           42 :                                                             la_min, lb_max, lb_min, npgfa, npgfb, &
     109           42 :                                                             nsgfa, nsgfb, row_blk_sizes
     110           42 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     111              :       INTEGER, DIMENSION(nexp_max)                       :: nct
     112              :       LOGICAL                                            :: found, new_atom_pair
     113              :       REAL(KIND=dp)                                      :: dab, dac, dbc, f0, radius
     114           42 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
     115           42 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: ppl_fwork, ppl_work
     116           42 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: hab, pab
     117              :       REAL(KIND=dp), DIMENSION(3)                        :: force_a, force_b, rab, rac, rbc
     118           42 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: alpha, set_radius_a, set_radius_b
     119           42 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ccval, cval, h_block, p_block, rpgfa, &
     120           42 :                                                             rpgfb, sphi_a, sphi_b, zeta, zetb
     121           42 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_kg
     122           42 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     123              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     124              :       TYPE(local_potential_type), POINTER                :: tnadd_potential
     125              :       TYPE(neighbor_list_iterator_p_type), &
     126           42 :          DIMENSION(:), POINTER                           :: ap_iterator, nl_iterator
     127              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     128           42 :          POINTER                                         :: sac_kin
     129              : 
     130           84 :       IF (calculate_forces) THEN
     131           20 :          CALL timeset(routineN//"_forces", handle)
     132              :       ELSE
     133           22 :          CALL timeset(routineN, handle)
     134              :       END IF
     135              : 
     136           42 :       NULLIFY (matrix_kg)
     137           42 :       IF (ASSOCIATED(kg_env%tnadd_mat)) THEN
     138           30 :          CALL dbcsr_deallocate_matrix_set(kg_env%tnadd_mat)
     139              :       END IF
     140           42 :       sac_kin => kg_env%sac_kin
     141           42 :       atom_to_molecule => kg_env%atom_to_molecule
     142              : 
     143           42 :       nkind = SIZE(atomic_kind_set)
     144           42 :       natom = SIZE(particle_set)
     145              : 
     146           42 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     147              : 
     148           42 :       IF (calculate_forces) THEN
     149           20 :          nder = 1
     150           20 :          IF (SIZE(matrix_p, 1) == 2) THEN
     151            0 :             DO img = 1, SIZE(matrix_p, 2)
     152              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     153            0 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     154              :                CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
     155            0 :                               alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
     156              :             END DO
     157              :          END IF
     158              :       ELSE
     159              :          nder = 0
     160              :       END IF
     161              : 
     162           42 :       maxder = ncoset(nder)
     163              : 
     164              :       CALL get_qs_kind_set(qs_kind_set, maxpol=maxpol, maxco=maxco, maxlgto=maxlgto, &
     165           42 :                            maxsgf=maxsgf, maxnset=maxnset)
     166              : 
     167           42 :       maxl = MAX(maxlgto, maxpol)
     168           42 :       CALL init_orbital_pointers(maxl + nder + 1)
     169              : 
     170           42 :       ldsab = MAX(maxco, ncoset(maxpol), maxsgf, maxpol)
     171           42 :       ldai = ncoset(maxl + nder + 1)
     172              : 
     173          192 :       ALLOCATE (basis_set_list(nkind))
     174          108 :       DO ikind = 1, nkind
     175           66 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
     176          108 :          IF (ASSOCIATED(basis_set_a)) THEN
     177           66 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     178              :          ELSE
     179            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     180              :          END IF
     181              :       END DO
     182              : 
     183              :       ! build the matrix
     184          168 :       ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
     185              : 
     186           42 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
     187           42 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes)
     188              : 
     189           42 :       CALL dbcsr_allocate_matrix_set(matrix_kg, 1)
     190              : 
     191           42 :       ALLOCATE (matrix_kg(1)%matrix)
     192              :       CALL dbcsr_create(matrix=matrix_kg(1)%matrix, &
     193              :                         name="Nonadditive kinetic energy potential", &
     194              :                         dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
     195              :                         row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, &
     196           42 :                         reuse_arrays=.TRUE.)
     197           42 :       CALL cp_dbcsr_alloc_block_from_nbl(matrix_kg(1)%matrix, sab_orb)
     198           42 :       CALL dbcsr_set(matrix_kg(1)%matrix, 0.0_dp)
     199              : 
     200              :       nthread = 1
     201           42 : !$    nthread = omp_get_max_threads()
     202              : 
     203           42 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb, nthread=nthread)
     204              :       ! iterator for basis/potential list
     205           42 :       CALL neighbor_list_iterator_create(ap_iterator, sac_kin, search=.TRUE., nthread=nthread)
     206              : 
     207              : !$OMP PARALLEL &
     208              : !$OMP DEFAULT (NONE) &
     209              : !$OMP SHARED  (nl_iterator, ap_iterator, basis_set_list, calculate_forces, force, use_virial,&
     210              : !$OMP          matrix_kg, matrix_p,virial, atomic_kind_set, qs_kind_set, particle_set,&
     211              : !$OMP          sab_orb, sac_kin, nthread, ncoset, nkind,&
     212              : !$OMP          atom_of_kind, ldsab,  maxnset, maxder, &
     213              : !$OMP          maxlgto, nder, maxco, atom_to_molecule) &
     214              : !$OMP PRIVATE (ikind, jkind, inode, iatom, jatom, rab, basis_set_a, basis_set_b, atom_a, &
     215              : !$OMP          atom_b, first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
     216              : !$OMP          zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
     217              : !$OMP          zetb, last_iatom, last_jatom, new_atom_pair, dab, irow, icol, h_block, found, iset, ncoa, &
     218              : !$OMP          sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
     219              : !$OMP          radius, alpha, nct, imol, jmol, kmol,&
     220              : !$OMP          npol, ngau, cval, ccval, rac, dac, rbc, dbc, &
     221              : !$OMP          set_radius_a,  rpgfa, force_a, force_b, ppl_fwork, mepos, &
     222              : !$OMP          f0, katom, ppl_work, atom_c,&
     223           42 : !$OMP          ldai,tnadd_potential)
     224              : 
     225              :       mepos = 0
     226              : !$    mepos = omp_get_thread_num()
     227              : 
     228              :       ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
     229              :       ldai = ncoset(2*maxlgto + 2*nder)
     230              :       ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2*nder + 1)))
     231              :       IF (calculate_forces) THEN
     232              :          ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
     233              :          ldai = ncoset(maxlgto)
     234              :          ALLOCATE (ppl_fwork(ldai, ldai, maxder))
     235              :       END IF
     236              : 
     237              :       last_iatom = 0
     238              :       last_jatom = 0
     239              :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     240              : 
     241              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, inode=inode, &
     242              :                                 iatom=iatom, jatom=jatom, r=rab)
     243              : 
     244              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     245              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     246              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     247              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     248              : 
     249              :          atom_a = atom_of_kind(iatom)
     250              :          atom_b = atom_of_kind(jatom)
     251              :          imol = atom_to_molecule(iatom)
     252              :          jmol = atom_to_molecule(jatom)
     253              :          CPASSERT(imol == jmol)
     254              : 
     255              :          ! basis ikind
     256              :          first_sgfa => basis_set_a%first_sgf
     257              :          la_max => basis_set_a%lmax
     258              :          la_min => basis_set_a%lmin
     259              :          npgfa => basis_set_a%npgf
     260              :          nseta = basis_set_a%nset
     261              :          nsgfa => basis_set_a%nsgf_set
     262              :          rpgfa => basis_set_a%pgf_radius
     263              :          set_radius_a => basis_set_a%set_radius
     264              :          sphi_a => basis_set_a%sphi
     265              :          zeta => basis_set_a%zet
     266              :          ! basis jkind
     267              :          first_sgfb => basis_set_b%first_sgf
     268              :          lb_max => basis_set_b%lmax
     269              :          lb_min => basis_set_b%lmin
     270              :          npgfb => basis_set_b%npgf
     271              :          nsetb = basis_set_b%nset
     272              :          nsgfb => basis_set_b%nsgf_set
     273              :          rpgfb => basis_set_b%pgf_radius
     274              :          set_radius_b => basis_set_b%set_radius
     275              :          sphi_b => basis_set_b%sphi
     276              :          zetb => basis_set_b%zet
     277              : 
     278              :          dab = SQRT(SUM(rab*rab))
     279              : 
     280              :          IF (iatom == last_iatom .AND. jatom == last_jatom) THEN
     281              :             new_atom_pair = .FALSE.
     282              :          ELSE
     283              :             new_atom_pair = .TRUE.
     284              :             last_jatom = jatom
     285              :             last_iatom = iatom
     286              :          END IF
     287              : 
     288              :          ! *** Use the symmetry of the first derivatives ***
     289              :          IF (iatom == jatom) THEN
     290              :             f0 = 1.0_dp
     291              :          ELSE
     292              :             f0 = 2.0_dp
     293              :          END IF
     294              : 
     295              :          ! *** Create matrix blocks for a new matrix block column ***
     296              :          IF (new_atom_pair) THEN
     297              :             IF (iatom <= jatom) THEN
     298              :                irow = iatom
     299              :                icol = jatom
     300              :             ELSE
     301              :                irow = jatom
     302              :                icol = iatom
     303              :             END IF
     304              :             NULLIFY (h_block)
     305              :             CALL dbcsr_get_block_p(matrix_kg(1)%matrix, irow, icol, h_block, found)
     306              :             IF (ASSOCIATED(h_block)) THEN
     307              :             IF (calculate_forces) THEN
     308              :                CPASSERT(SIZE(matrix_p, 2) == 1)
     309              :                NULLIFY (p_block)
     310              :                CALL dbcsr_get_block_p(matrix_p(1, 1)%matrix, irow, icol, p_block, found)
     311              :                IF (ASSOCIATED(p_block)) THEN
     312              :                   DO iset = 1, nseta
     313              :                      ncoa = npgfa(iset)*ncoset(la_max(iset))
     314              :                      sgfa = first_sgfa(1, iset)
     315              :                      DO jset = 1, nsetb
     316              :                         ncob = npgfb(jset)*ncoset(lb_max(jset))
     317              :                         sgfb = first_sgfb(1, jset)
     318              : 
     319              :                         ! *** Decontract density matrix block ***
     320              :                         IF (iatom <= jatom) THEN
     321              :                            CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
     322              :                                       1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     323              :                                       p_block(sgfa, sgfb), SIZE(p_block, 1), &
     324              :                                       0.0_dp, work(1, 1), SIZE(work, 1))
     325              :                         ELSE
     326              :                            CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), &
     327              :                                       1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     328              :                                       p_block(sgfb, sgfa), SIZE(p_block, 1), &
     329              :                                       0.0_dp, work(1, 1), SIZE(work, 1))
     330              :                         END IF
     331              : 
     332              :                         CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
     333              :                                    1.0_dp, work(1, 1), SIZE(work, 1), &
     334              :                                    sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     335              :                                    0.0_dp, pab(1, 1, iset, jset), SIZE(pab, 1))
     336              :                      END DO
     337              :                   END DO
     338              :                END IF
     339              :             END IF
     340              :             END IF
     341              :          END IF
     342              : 
     343              :          hab = 0._dp
     344              : 
     345              :          ! loop over all kinds for nonadditive kinetic energy potential atoms
     346              :          DO kkind = 1, nkind
     347              : 
     348              :             CALL get_qs_kind(qs_kind_set(kkind), tnadd_potential=tnadd_potential)
     349              :             IF (.NOT. ASSOCIATED(tnadd_potential)) CYCLE
     350              :             CALL get_potential(potential=tnadd_potential, &
     351              :                                alpha=alpha, cval=cval, ngau=ngau, npol=npol, radius=radius)
     352              :             nct = npol
     353              :             ALLOCATE (ccval(npol, ngau))
     354              :             ccval(1:npol, 1:ngau) = TRANSPOSE(cval(1:ngau, 1:npol))
     355              : 
     356              :             CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos)
     357              : 
     358              :             DO WHILE (nl_sub_iterate(ap_iterator, mepos) == 0)
     359              : 
     360              :                CALL get_iterator_info(ap_iterator, mepos, jatom=katom, r=rac)
     361              : 
     362              :                ! this potential only acts on other moleclules
     363              :                kmol = atom_to_molecule(katom)
     364              :                IF (kmol == imol) CYCLE
     365              : 
     366              :                dac = SQRT(SUM(rac*rac))
     367              :                rbc(:) = rac(:) - rab(:)
     368              :                dbc = SQRT(SUM(rbc*rbc))
     369              :                IF ((MAXVAL(set_radius_a(:)) + radius < dac) .OR. &
     370              :                    (MAXVAL(set_radius_b(:)) + radius < dbc)) THEN
     371              :                   CYCLE
     372              :                END IF
     373              : 
     374              :                DO iset = 1, nseta
     375              :                   IF (set_radius_a(iset) + radius < dac) CYCLE
     376              :                   ncoa = npgfa(iset)*ncoset(la_max(iset))
     377              :                   sgfa = first_sgfa(1, iset)
     378              :                   DO jset = 1, nsetb
     379              :                      IF (set_radius_b(jset) + radius < dbc) CYCLE
     380              :                      ncob = npgfb(jset)*ncoset(lb_max(jset))
     381              :                      sgfb = first_sgfb(1, jset)
     382              :                      IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     383              :                      ! *** Calculate the GTH pseudo potential forces ***
     384              :                      IF (calculate_forces) THEN
     385              : 
     386              :                         CALL ppl_integral( &
     387              :                            la_max(iset), la_min(iset), npgfa(iset), &
     388              :                            rpgfa(:, iset), zeta(:, iset), &
     389              :                            lb_max(jset), lb_min(jset), npgfb(jset), &
     390              :                            rpgfb(:, jset), zetb(:, jset), &
     391              :                            ngau, alpha, nct, ccval, radius, &
     392              :                            rab, dab, rac, dac, rbc, dbc, &
     393              :                            hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
     394              :                            force_a, force_b, ppl_fwork)
     395              :                         ! *** The derivatives w.r.t. atomic center c are    ***
     396              :                         ! *** calculated using the translational invariance ***
     397              :                         ! *** of the first derivatives                      ***
     398              :                         atom_c = atom_of_kind(katom)
     399              : 
     400              : !$OMP CRITICAL(force_critical)
     401              :                         force(ikind)%kinetic(1, atom_a) = force(ikind)%kinetic(1, atom_a) + f0*force_a(1)
     402              :                         force(ikind)%kinetic(2, atom_a) = force(ikind)%kinetic(2, atom_a) + f0*force_a(2)
     403              :                         force(ikind)%kinetic(3, atom_a) = force(ikind)%kinetic(3, atom_a) + f0*force_a(3)
     404              :                         force(kkind)%kinetic(1, atom_c) = force(kkind)%kinetic(1, atom_c) - f0*force_a(1)
     405              :                         force(kkind)%kinetic(2, atom_c) = force(kkind)%kinetic(2, atom_c) - f0*force_a(2)
     406              :                         force(kkind)%kinetic(3, atom_c) = force(kkind)%kinetic(3, atom_c) - f0*force_a(3)
     407              : 
     408              :                         force(jkind)%kinetic(1, atom_b) = force(jkind)%kinetic(1, atom_b) + f0*force_b(1)
     409              :                         force(jkind)%kinetic(2, atom_b) = force(jkind)%kinetic(2, atom_b) + f0*force_b(2)
     410              :                         force(jkind)%kinetic(3, atom_b) = force(jkind)%kinetic(3, atom_b) + f0*force_b(3)
     411              :                         force(kkind)%kinetic(1, atom_c) = force(kkind)%kinetic(1, atom_c) - f0*force_b(1)
     412              :                         force(kkind)%kinetic(2, atom_c) = force(kkind)%kinetic(2, atom_c) - f0*force_b(2)
     413              :                         force(kkind)%kinetic(3, atom_c) = force(kkind)%kinetic(3, atom_c) - f0*force_b(3)
     414              : 
     415              :                         IF (use_virial) THEN
     416              :                            CALL virial_pair_force(virial%pv_virial, f0, force_a, rac)
     417              :                            CALL virial_pair_force(virial%pv_virial, f0, force_b, rbc)
     418              :                         END IF
     419              : !$OMP END CRITICAL(force_critical)
     420              : 
     421              :                      ELSE
     422              :                         CALL ppl_integral( &
     423              :                            la_max(iset), la_min(iset), npgfa(iset), &
     424              :                            rpgfa(:, iset), zeta(:, iset), &
     425              :                            lb_max(jset), lb_min(jset), npgfb(jset), &
     426              :                            rpgfb(:, jset), zetb(:, jset), &
     427              :                            ngau, alpha, nct, ccval, radius, &
     428              :                            rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
     429              :                      END IF
     430              :                   END DO
     431              :                END DO
     432              :             END DO
     433              :             DEALLOCATE (ccval)
     434              :          END DO
     435              : 
     436              :          ! *** Contract integrals
     437              :          DO iset = 1, nseta
     438              :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     439              :             sgfa = first_sgfa(1, iset)
     440              :             DO jset = 1, nsetb
     441              :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     442              :                sgfb = first_sgfb(1, jset)
     443              : 
     444              :                CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     445              :                           1.0_dp, hab(1, 1, iset, jset), SIZE(hab, 1), &
     446              :                           sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     447              :                           0.0_dp, work(1, 1), SIZE(work, 1))
     448              : 
     449              : !$OMP CRITICAL(h_block_critical)
     450              :                IF (iatom <= jatom) THEN
     451              :                   CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     452              :                              1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     453              :                              work(1, 1), SIZE(work, 1), &
     454              :                              1.0_dp, h_block(sgfa, sgfb), SIZE(h_block, 1))
     455              :                ELSE
     456              :                   CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
     457              :                              1.0_dp, work(1, 1), SIZE(work, 1), &
     458              :                              sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     459              :                              1.0_dp, h_block(sgfb, sgfa), SIZE(h_block, 1))
     460              :                END IF
     461              : !$OMP END CRITICAL(h_block_critical)
     462              : 
     463              :             END DO
     464              :          END DO
     465              :       END DO
     466              : 
     467              :       DEALLOCATE (hab, work, ppl_work)
     468              : 
     469              :       IF (calculate_forces) THEN
     470              :          DEALLOCATE (pab, ppl_fwork)
     471              :       END IF
     472              : 
     473              : !$OMP END PARALLEL
     474              : 
     475           42 :       CALL neighbor_list_iterator_release(ap_iterator)
     476           42 :       CALL neighbor_list_iterator_release(nl_iterator)
     477              : 
     478           84 :       DO i = 1, SIZE(matrix_kg)
     479           84 :          CALL dbcsr_finalize(matrix_kg(i)%matrix)
     480              :       END DO
     481           42 :       kg_env%tnadd_mat => matrix_kg
     482              : 
     483           42 :       DEALLOCATE (basis_set_list)
     484              : 
     485           42 :       IF (calculate_forces) THEN
     486              :          ! *** If LSD, then recover alpha density and beta density     ***
     487              :          ! *** from the total density (1) and the spin density (2)     ***
     488           20 :          IF (SIZE(matrix_p, 1) == 2) THEN
     489            0 :             DO img = 1, SIZE(matrix_p, 2)
     490              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     491            0 :                               alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
     492              :                CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
     493            0 :                               alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
     494              :             END DO
     495              :          END IF
     496              :       END IF
     497              : 
     498           42 :       CALL timestop(handle)
     499              : 
     500          126 :    END SUBROUTINE build_tnadd_mat
     501              : 
     502              : !==========================================================================================================
     503              : 
     504              : END MODULE kg_tnadd_mat
        

Generated by: LCOV version 2.0-1