LCOV - code coverage report
Current view: top level - src - qs_elec_field.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.1 % 113 112
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              : ! **************************************************************************************************
       9              : !> \brief Distribution of the electric field gradient integral matrix.
      10              : !> \par History
      11              : !> \author VW (27.02.2009)
      12              : ! **************************************************************************************************
      13              : MODULE qs_elec_field
      14              : 
      15              :    USE ai_elec_field,                   ONLY: efg
      16              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      17              :                                               gto_basis_set_type
      18              :    USE block_p_types,                   ONLY: block_p_type
      19              :    USE cell_types,                      ONLY: cell_type,&
      20              :                                               pbc
      21              :    USE cp_dbcsr_api,                    ONLY: dbcsr_get_block_p,&
      22              :                                               dbcsr_p_type
      23              :    USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_sparse_matrix
      24              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      25              :                                               cp_logger_type
      26              :    USE cp_output_handling,              ONLY: cp_p_file,&
      27              :                                               cp_print_key_finished_output,&
      28              :                                               cp_print_key_should_output,&
      29              :                                               cp_print_key_unit_nr
      30              :    USE input_section_types,             ONLY: section_vals_val_get
      31              :    USE kinds,                           ONLY: dp
      32              :    USE message_passing,                 ONLY: mp_para_env_type
      33              :    USE orbital_pointers,                ONLY: init_orbital_pointers,&
      34              :                                               ncoset
      35              :    USE particle_types,                  ONLY: particle_type
      36              :    USE qs_environment_types,            ONLY: get_qs_env,&
      37              :                                               qs_environment_type
      38              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      39              :                                               get_qs_kind_set,&
      40              :                                               qs_kind_type
      41              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      42              :                                               neighbor_list_iterate,&
      43              :                                               neighbor_list_iterator_create,&
      44              :                                               neighbor_list_iterator_p_type,&
      45              :                                               neighbor_list_iterator_release,&
      46              :                                               neighbor_list_set_p_type
      47              : #include "./base/base_uses.f90"
      48              : 
      49              :    IMPLICIT NONE
      50              : 
      51              :    PRIVATE
      52              : 
      53              : ! *** Global parameters ***
      54              : 
      55              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_elec_field'
      56              : 
      57              : ! *** Public subroutines ***
      58              : 
      59              :    PUBLIC :: build_efg_matrix
      60              : 
      61              : CONTAINS
      62              : 
      63              : ! **************************************************************************************************
      64              : !> \brief   Calculation of the electric field gradient matrix over
      65              : !>          Cartesian Gaussian functions.
      66              : !> \param qs_env ...
      67              : !> \param matrix_efg ...
      68              : !> \param rc ...
      69              : !> \date    27.02.2009
      70              : !> \author  VW
      71              : !> \version 1.0
      72              : ! **************************************************************************************************
      73              : 
      74           44 :    SUBROUTINE build_efg_matrix(qs_env, matrix_efg, rc)
      75              : 
      76              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      77              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_efg
      78              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rc
      79              : 
      80              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'build_efg_matrix'
      81              : 
      82              :       INTEGER :: after, handle, i, iatom, icol, ikind, inode, irow, iset, iw, jatom, jkind, jset, &
      83              :          last_jatom, ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, &
      84              :          nseta, nsetb, sgfa, sgfb
      85           44 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
      86           44 :                                                             npgfb, nsgfa, nsgfb
      87           44 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
      88              :       LOGICAL                                            :: found, new_atom_b, omit_headers
      89              :       REAL(KIND=dp)                                      :: dab, rab2
      90           44 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: work
      91           44 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: efgab, rr_work
      92              :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rab, rac, rb, rbc
      93           44 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
      94           44 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
      95           44 :       TYPE(block_p_type), ALLOCATABLE, DIMENSION(:)      :: efgint
      96              :       TYPE(cell_type), POINTER                           :: cell
      97              :       TYPE(cp_logger_type), POINTER                      :: logger
      98           44 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
      99              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     100              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     101              :       TYPE(neighbor_list_iterator_p_type), &
     102           44 :          DIMENSION(:), POINTER                           :: nl_iterator
     103              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     104           44 :          POINTER                                         :: sab_orb
     105           44 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     106           44 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     107              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     108              : 
     109           44 :       CALL timeset(routineN, handle)
     110              : 
     111           44 :       NULLIFY (cell, sab_orb, qs_kind_set, particle_set, para_env)
     112           44 :       NULLIFY (logger)
     113              : 
     114           44 :       logger => cp_get_default_logger()
     115              : 
     116              :       CALL get_qs_env(qs_env=qs_env, &
     117              :                       qs_kind_set=qs_kind_set, &
     118              :                       particle_set=particle_set, &
     119              :                       neighbor_list_id=neighbor_list_id, &
     120              :                       para_env=para_env, &
     121              :                       sab_orb=sab_orb, &
     122           44 :                       cell=cell)
     123              : 
     124           44 :       nkind = SIZE(qs_kind_set)
     125           44 :       natom = SIZE(particle_set)
     126              : 
     127              :       ! *** Allocate work storage ***
     128              :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     129              :                            maxco=maxco, &
     130              :                            maxlgto=maxlgto, &
     131           44 :                            maxsgf=maxsgf)
     132              : 
     133           44 :       ldai = ncoset(maxlgto + 2)
     134           44 :       CALL init_orbital_pointers(ldai)
     135              : 
     136          220 :       ALLOCATE (rr_work(0:2*maxlgto + 4, ldai, ldai))
     137              : 
     138          220 :       ALLOCATE (efgab(maxco, maxco, 6))
     139              : 
     140          176 :       ALLOCATE (work(maxco, maxsgf))
     141              : 
     142          308 :       ALLOCATE (efgint(6))
     143              : 
     144       141724 :       rr_work(:, :, :) = 0.0_dp
     145        18692 :       efgab(:, :, :) = 0.0_dp
     146         3496 :       work(:, :) = 0.0_dp
     147              : 
     148          218 :       ALLOCATE (basis_set_list(nkind))
     149          130 :       DO ikind = 1, nkind
     150           86 :          qs_kind => qs_kind_set(ikind)
     151           86 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
     152          130 :          IF (ASSOCIATED(basis_set_a)) THEN
     153           86 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     154              :          ELSE
     155            0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     156              :          END IF
     157              :       END DO
     158           44 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     159         4918 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     160              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
     161         4874 :                                 iatom=iatom, jatom=jatom, r=rab)
     162         4874 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     163         4874 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     164         4874 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     165         4874 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     166         4874 :          ra = pbc(particle_set(iatom)%r, cell)
     167              :          ! basis ikind
     168         4874 :          first_sgfa => basis_set_a%first_sgf
     169         4874 :          la_max => basis_set_a%lmax
     170         4874 :          la_min => basis_set_a%lmin
     171         4874 :          npgfa => basis_set_a%npgf
     172         4874 :          nseta = basis_set_a%nset
     173         4874 :          nsgfa => basis_set_a%nsgf_set
     174         4874 :          rpgfa => basis_set_a%pgf_radius
     175         4874 :          set_radius_a => basis_set_a%set_radius
     176         4874 :          sphi_a => basis_set_a%sphi
     177         4874 :          zeta => basis_set_a%zet
     178              :          ! basis jkind
     179         4874 :          first_sgfb => basis_set_b%first_sgf
     180         4874 :          lb_max => basis_set_b%lmax
     181         4874 :          lb_min => basis_set_b%lmin
     182         4874 :          npgfb => basis_set_b%npgf
     183         4874 :          nsetb = basis_set_b%nset
     184         4874 :          nsgfb => basis_set_b%nsgf_set
     185         4874 :          rpgfb => basis_set_b%pgf_radius
     186         4874 :          set_radius_b => basis_set_b%set_radius
     187         4874 :          sphi_b => basis_set_b%sphi
     188         4874 :          zetb => basis_set_b%zet
     189              : 
     190         4874 :          IF (inode == 1) last_jatom = 0
     191              : 
     192        19496 :          rb = rab + ra
     193         4874 :          rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     194         4874 :          dab = SQRT(rab2)
     195         4874 :          rac = pbc(ra, rc, cell)
     196        19496 :          rbc = rac - rab
     197              : 
     198         4874 :          IF (jatom /= last_jatom) THEN
     199              :             new_atom_b = .TRUE.
     200              :             last_jatom = jatom
     201              :          ELSE
     202              :             new_atom_b = .FALSE.
     203              :          END IF
     204              : 
     205              :          IF (new_atom_b) THEN
     206          272 :             IF (iatom <= jatom) THEN
     207          164 :                irow = iatom
     208          164 :                icol = jatom
     209              :             ELSE
     210          108 :                irow = jatom
     211          108 :                icol = iatom
     212              :             END IF
     213              : 
     214         1904 :             DO i = 1, 6
     215         1632 :                NULLIFY (efgint(i)%block)
     216              :                CALL dbcsr_get_block_p(matrix=matrix_efg(i)%matrix, &
     217         1904 :                                       row=irow, col=icol, BLOCK=efgint(i)%block, found=found)
     218              :             END DO
     219              :          END IF
     220              : 
     221        15700 :          DO iset = 1, nseta
     222              : 
     223        10782 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     224        10782 :             sgfa = first_sgfa(1, iset)
     225              : 
     226        39798 :             DO jset = 1, nsetb
     227              : 
     228        24142 :                IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     229              : 
     230         9222 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     231         9222 :                sgfb = first_sgfb(1, jset)
     232              : 
     233              :                ! *** Calculate the primitive fermi contact integrals ***
     234              : 
     235              :                CALL efg(la_max(iset), la_min(iset), npgfa(iset), &
     236              :                         rpgfa(:, iset), zeta(:, iset), &
     237              :                         lb_max(jset), lb_min(jset), npgfb(jset), &
     238              :                         rpgfb(:, jset), zetb(:, jset), &
     239         9222 :                         rac, rbc, rab, efgab, SIZE(rr_work, 1), SIZE(rr_work, 2), rr_work)
     240              : 
     241              :                ! *** Contraction step ***
     242              : 
     243        75336 :                DO i = 1, 6
     244              : 
     245              :                   CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
     246              :                              1.0_dp, efgab(1, 1, i), SIZE(efgab, 1), &
     247              :                              sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     248        55332 :                              0.0_dp, work(1, 1), SIZE(work, 1))
     249              : 
     250        79474 :                   IF (iatom <= jatom) THEN
     251              :                      CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
     252              :                                 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     253              :                                 work(1, 1), SIZE(work, 1), &
     254              :                                 1.0_dp, efgint(i)%block(sgfa, sgfb), &
     255        33582 :                                 SIZE(efgint(i)%block, 1))
     256              : 
     257              :                   ELSE
     258              : 
     259              :                      CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
     260              :                                 1.0_dp, work(1, 1), SIZE(work, 1), &
     261              :                                 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     262              :                                 1.0_dp, efgint(i)%block(sgfb, sgfa), &
     263        21750 :                                 SIZE(efgint(i)%block, 1))
     264              :                   END IF
     265              : 
     266              :                END DO
     267              : 
     268              :             END DO
     269              : 
     270              :          END DO
     271              : 
     272              :       END DO
     273           44 :       CALL neighbor_list_iterator_release(nl_iterator)
     274              : 
     275              :       ! *** Release work storage ***
     276           44 :       DEALLOCATE (basis_set_list)
     277              : 
     278           44 :       DEALLOCATE (efgab)
     279              : 
     280           44 :       DEALLOCATE (work)
     281              : 
     282           44 :       DEALLOCATE (efgint)
     283              : 
     284              :       ! Print the electric field gradient matrix, if requested
     285              : 
     286           44 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     287              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/EFG"), cp_p_file)) THEN
     288              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/EFG", &
     289            8 :                                    extension=".Log")
     290            8 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
     291            8 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
     292            8 :          after = MIN(MAX(after, 1), 16)
     293              :          CALL cp_dbcsr_write_sparse_matrix(matrix_efg(1)%matrix, 4, after, qs_env, &
     294            8 :                                            para_env, output_unit=iw, omit_headers=omit_headers)
     295              :          CALL cp_dbcsr_write_sparse_matrix(matrix_efg(2)%matrix, 4, after, qs_env, &
     296            8 :                                            para_env, output_unit=iw, omit_headers=omit_headers)
     297              :          CALL cp_dbcsr_write_sparse_matrix(matrix_efg(3)%matrix, 4, after, qs_env, &
     298            8 :                                            para_env, output_unit=iw, omit_headers=omit_headers)
     299              :          CALL cp_dbcsr_write_sparse_matrix(matrix_efg(4)%matrix, 4, after, qs_env, &
     300            8 :                                            para_env, output_unit=iw, omit_headers=omit_headers)
     301              :          CALL cp_dbcsr_write_sparse_matrix(matrix_efg(5)%matrix, 4, after, qs_env, &
     302            8 :                                            para_env, output_unit=iw, omit_headers=omit_headers)
     303              :          CALL cp_dbcsr_write_sparse_matrix(matrix_efg(6)%matrix, 4, after, qs_env, &
     304            8 :                                            para_env, output_unit=iw, omit_headers=omit_headers)
     305              :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
     306            8 :                                            "DFT%PRINT%AO_MATRICES/EFG")
     307              :       END IF
     308              : 
     309           44 :       CALL timestop(handle)
     310              : 
     311          132 :    END SUBROUTINE build_efg_matrix
     312              : 
     313              : ! **************************************************************************************************
     314              : 
     315              : END MODULE qs_elec_field
     316              : 
        

Generated by: LCOV version 2.0-1