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

            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 Calculates the energy contribution and the mo_derivative of
      10              : !>        a static electric field (nonperiodic)
      11              : !> \par History
      12              : !>      Adjusted from qs_efield_local
      13              : !> \author JGH (10.2019)
      14              : ! **************************************************************************************************
      15              : MODULE ec_efield_local
      16              :    USE ai_moments,                      ONLY: dipole_force
      17              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18              :                                               get_atomic_kind,&
      19              :                                               get_atomic_kind_set
      20              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      21              :                                               gto_basis_set_type
      22              :    USE cell_types,                      ONLY: cell_type,&
      23              :                                               pbc
      24              :    USE cp_control_types,                ONLY: dft_control_type
      25              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      26              :                                               dbcsr_copy,&
      27              :                                               dbcsr_get_block_p,&
      28              :                                               dbcsr_p_type,&
      29              :                                               dbcsr_set
      30              :    USE ec_env_types,                    ONLY: energy_correction_type
      31              :    USE kinds,                           ONLY: dp
      32              :    USE message_passing,                 ONLY: mp_para_env_type
      33              :    USE orbital_pointers,                ONLY: ncoset
      34              :    USE particle_types,                  ONLY: particle_type
      35              :    USE qs_energy_types,                 ONLY: qs_energy_type
      36              :    USE qs_environment_types,            ONLY: get_qs_env,&
      37              :                                               qs_environment_type
      38              :    USE qs_force_types,                  ONLY: qs_force_type
      39              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      40              :                                               qs_kind_type
      41              :    USE qs_moments,                      ONLY: build_local_moment_matrix
      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              :    USE qs_period_efield_types,          ONLY: efield_berry_type,&
      49              :                                               init_efield_matrices,&
      50              :                                               set_efield_matrices
      51              : #include "./base/base_uses.f90"
      52              : 
      53              :    IMPLICIT NONE
      54              : 
      55              :    PRIVATE
      56              : 
      57              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_efield_local'
      58              : 
      59              :    ! *** Public subroutines ***
      60              : 
      61              :    PUBLIC :: ec_efield_local_operator, ec_efield_integrals
      62              : 
      63              : ! **************************************************************************************************
      64              : 
      65              : CONTAINS
      66              : 
      67              : ! **************************************************************************************************
      68              : 
      69              : ! **************************************************************************************************
      70              : !> \brief ...
      71              : !> \param qs_env ...
      72              : !> \param ec_env ...
      73              : !> \param calculate_forces ...
      74              : ! **************************************************************************************************
      75         1126 :    SUBROUTINE ec_efield_local_operator(qs_env, ec_env, calculate_forces)
      76              : 
      77              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      78              :       TYPE(energy_correction_type), POINTER              :: ec_env
      79              :       LOGICAL, INTENT(IN)                                :: calculate_forces
      80              : 
      81              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_local_operator'
      82              : 
      83              :       INTEGER                                            :: handle
      84              :       REAL(dp), DIMENSION(3)                             :: rpoint
      85              :       TYPE(dft_control_type), POINTER                    :: dft_control
      86              : 
      87         1126 :       CALL timeset(routineN, handle)
      88              : 
      89         1126 :       NULLIFY (dft_control)
      90         1126 :       CALL get_qs_env(qs_env, dft_control=dft_control)
      91              : 
      92         1126 :       IF (dft_control%apply_efield) THEN
      93           86 :          rpoint = 0.0_dp
      94           86 :          CALL ec_efield_integrals(qs_env, ec_env, rpoint)
      95           86 :          CALL ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces)
      96              :       END IF
      97              : 
      98         1126 :       CALL timestop(handle)
      99              : 
     100         1126 :    END SUBROUTINE ec_efield_local_operator
     101              : 
     102              : ! **************************************************************************************************
     103              : !> \brief ...
     104              : !> \param qs_env ...
     105              : !> \param ec_env ...
     106              : !> \param rpoint ...
     107              : ! **************************************************************************************************
     108          106 :    SUBROUTINE ec_efield_integrals(qs_env, ec_env, rpoint)
     109              : 
     110              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     111              :       TYPE(energy_correction_type), POINTER              :: ec_env
     112              :       REAL(dp), DIMENSION(3), INTENT(IN)                 :: rpoint
     113              : 
     114              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_integrals'
     115              : 
     116              :       INTEGER                                            :: handle, i
     117          106 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: dipmat, matrix_s
     118              :       TYPE(efield_berry_type), POINTER                   :: efield, efieldref
     119              : 
     120          106 :       CALL timeset(routineN, handle)
     121              : 
     122          106 :       CALL get_qs_env(qs_env=qs_env, efield=efieldref)
     123          106 :       efield => ec_env%efield
     124          106 :       CALL init_efield_matrices(efield)
     125          106 :       matrix_s => ec_env%matrix_s(:, 1)
     126          424 :       ALLOCATE (dipmat(3))
     127          424 :       DO i = 1, 3
     128          318 :          ALLOCATE (dipmat(i)%matrix)
     129          318 :          CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'DIP MAT')
     130          424 :          CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
     131              :       END DO
     132          106 :       CALL build_local_moment_matrix(qs_env, dipmat, 1, rpoint, basis_type="HARRIS")
     133          106 :       CALL set_efield_matrices(efield=efield, dipmat=dipmat)
     134          106 :       ec_env%efield => efield
     135              : 
     136          106 :       CALL timestop(handle)
     137              : 
     138          106 :    END SUBROUTINE ec_efield_integrals
     139              : 
     140              : ! **************************************************************************************************
     141              : !> \brief ...
     142              : !> \param qs_env ...
     143              : !> \param ec_env ...
     144              : !> \param rpoint ...
     145              : !> \param calculate_forces ...
     146              : ! **************************************************************************************************
     147           86 :    SUBROUTINE ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces)
     148              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     149              :       TYPE(energy_correction_type), POINTER              :: ec_env
     150              :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rpoint
     151              :       LOGICAL                                            :: calculate_forces
     152              : 
     153              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_efield_mo_derivatives'
     154              : 
     155              :       INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, &
     156              :          jatom, jkind, jset, ldab, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
     157           86 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     158           86 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     159           86 :                                                             npgfb, nsgfa, nsgfb
     160           86 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     161              :       LOGICAL                                            :: found, trans
     162              :       REAL(dp)                                           :: charge, dab, fdir
     163              :       REAL(dp), DIMENSION(3)                             :: ci, fieldpol, ra, rab, rac, rbc, ria
     164              :       REAL(dp), DIMENSION(3, 3)                          :: forcea, forceb
     165           86 :       REAL(dp), DIMENSION(:, :), POINTER                 :: p_block_a, p_block_b, pblock, pmat, work
     166           86 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     167           86 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
     168           86 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     169              :       TYPE(cell_type), POINTER                           :: cell
     170           86 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: dipmat, matrix_ks
     171              :       TYPE(dft_control_type), POINTER                    :: dft_control
     172              :       TYPE(efield_berry_type), POINTER                   :: efield
     173           86 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     174              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     175              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     176              :       TYPE(neighbor_list_iterator_p_type), &
     177           86 :          DIMENSION(:), POINTER                           :: nl_iterator
     178              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     179           86 :          POINTER                                         :: sab_orb
     180           86 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     181              :       TYPE(qs_energy_type), POINTER                      :: energy
     182           86 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     183           86 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     184              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     185              : 
     186           86 :       CALL timeset(routineN, handle)
     187              : 
     188           86 :       CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
     189              :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
     190           86 :                       energy=energy, para_env=para_env, sab_orb=sab_orb)
     191              : 
     192           86 :       efield => ec_env%efield
     193              : 
     194              :       fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
     195          344 :                  dft_control%efield_fields(1)%efield%strength
     196              : 
     197              :       ! nuclear contribution
     198           86 :       natom = SIZE(particle_set)
     199           86 :       IF (calculate_forces) THEN
     200           14 :          CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
     201           14 :          CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
     202              :       END IF
     203           86 :       ci = 0.0_dp
     204          340 :       DO ia = 1, natom
     205          254 :          CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
     206          254 :          CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
     207         1016 :          ria = particle_set(ia)%r - rpoint
     208         1016 :          ria = pbc(ria, cell)
     209         1016 :          ci(:) = ci(:) + charge*ria(:)
     210          594 :          IF (calculate_forces) THEN
     211           40 :             IF (para_env%mepos == 0) THEN
     212           20 :                iatom = atom_of_kind(ia)
     213           80 :                DO idir = 1, 3
     214           80 :                   force(ikind)%efield(idir, iatom) = force(ikind)%efield(idir, iatom) - fieldpol(idir)*charge
     215              :                END DO
     216              :             END IF
     217              :          END IF
     218              :       END DO
     219              : 
     220           86 :       IF (ec_env%should_update) THEN
     221          280 :          ec_env%efield_nuclear = -SUM(ci(:)*fieldpol(:))
     222              :          ! Update KS matrix
     223           70 :          matrix_ks => ec_env%matrix_h(:, 1)
     224           70 :          dipmat => efield%dipmat
     225          140 :          DO ispin = 1, SIZE(matrix_ks)
     226          350 :             DO idir = 1, 3
     227              :                CALL dbcsr_add(matrix_ks(ispin)%matrix, dipmat(idir)%matrix, &
     228          280 :                               alpha_scalar=1.0_dp, beta_scalar=fieldpol(idir))
     229              :             END DO
     230              :          END DO
     231              :       END IF
     232              : 
     233              :       ! forces from the efield contribution
     234           86 :       IF (calculate_forces) THEN
     235           14 :          nkind = SIZE(qs_kind_set)
     236           14 :          natom = SIZE(particle_set)
     237              : 
     238           70 :          ALLOCATE (basis_set_list(nkind))
     239           42 :          DO ikind = 1, nkind
     240           28 :             qs_kind => qs_kind_set(ikind)
     241           28 :             CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type="HARRIS")
     242           42 :             IF (ASSOCIATED(basis_set_a)) THEN
     243           28 :                basis_set_list(ikind)%gto_basis_set => basis_set_a
     244              :             ELSE
     245            0 :                NULLIFY (basis_set_list(ikind)%gto_basis_set)
     246              :             END IF
     247              :          END DO
     248              :          !
     249           14 :          CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     250           91 :          DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     251              :             CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     252           77 :                                    iatom=iatom, jatom=jatom, r=rab)
     253           77 :             basis_set_a => basis_set_list(ikind)%gto_basis_set
     254           77 :             IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     255           77 :             basis_set_b => basis_set_list(jkind)%gto_basis_set
     256           77 :             IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     257              :             ! basis ikind
     258           77 :             first_sgfa => basis_set_a%first_sgf
     259           77 :             la_max => basis_set_a%lmax
     260           77 :             la_min => basis_set_a%lmin
     261           77 :             npgfa => basis_set_a%npgf
     262           77 :             nseta = basis_set_a%nset
     263           77 :             nsgfa => basis_set_a%nsgf_set
     264           77 :             rpgfa => basis_set_a%pgf_radius
     265           77 :             set_radius_a => basis_set_a%set_radius
     266           77 :             sphi_a => basis_set_a%sphi
     267           77 :             zeta => basis_set_a%zet
     268              :             ! basis jkind
     269           77 :             first_sgfb => basis_set_b%first_sgf
     270           77 :             lb_max => basis_set_b%lmax
     271           77 :             lb_min => basis_set_b%lmin
     272           77 :             npgfb => basis_set_b%npgf
     273           77 :             nsetb = basis_set_b%nset
     274           77 :             nsgfb => basis_set_b%nsgf_set
     275           77 :             rpgfb => basis_set_b%pgf_radius
     276           77 :             set_radius_b => basis_set_b%set_radius
     277           77 :             sphi_b => basis_set_b%sphi
     278           77 :             zetb => basis_set_b%zet
     279              : 
     280           77 :             atom_a = atom_of_kind(iatom)
     281           77 :             atom_b = atom_of_kind(jatom)
     282              : 
     283          308 :             ra(:) = particle_set(iatom)%r(:) - rpoint(:)
     284           77 :             rac(:) = pbc(ra(:), cell)
     285          308 :             rbc(:) = rac(:) + rab(:)
     286           77 :             dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
     287              : 
     288           77 :             IF (iatom <= jatom) THEN
     289           50 :                irow = iatom
     290           50 :                icol = jatom
     291           50 :                trans = .FALSE.
     292              :             ELSE
     293           27 :                irow = jatom
     294           27 :                icol = iatom
     295           27 :                trans = .TRUE.
     296              :             END IF
     297              : 
     298           77 :             fdir = 2.0_dp
     299           77 :             IF (iatom == jatom .AND. dab < 1.e-10_dp) fdir = 1.0_dp
     300              : 
     301              :             ! density matrix
     302           77 :             NULLIFY (p_block_a)
     303           77 :             CALL dbcsr_get_block_p(ec_env%matrix_p(1, 1)%matrix, irow, icol, p_block_a, found)
     304           77 :             IF (.NOT. found) CYCLE
     305           77 :             IF (SIZE(ec_env%matrix_p, 1) > 1) THEN
     306            0 :                NULLIFY (p_block_b)
     307            0 :                CALL dbcsr_get_block_p(ec_env%matrix_p(2, 1)%matrix, irow, icol, p_block_b, found)
     308            0 :                CPASSERT(found)
     309              :             END IF
     310           77 :             forcea = 0.0_dp
     311           77 :             forceb = 0.0_dp
     312              : 
     313          231 :             DO iset = 1, nseta
     314          154 :                ncoa = npgfa(iset)*ncoset(la_max(iset))
     315          154 :                sgfa = first_sgfa(1, iset)
     316          539 :                DO jset = 1, nsetb
     317          308 :                   IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     318          230 :                   ncob = npgfb(jset)*ncoset(lb_max(jset))
     319          230 :                   sgfb = first_sgfb(1, jset)
     320              :                   ! Calculate the primitive integrals (da|O|b) and (a|O|db)
     321          230 :                   ldab = MAX(ncoa, ncob)
     322         1610 :                   ALLOCATE (work(ldab, ldab), pmat(ncoa, ncob))
     323              :                   ! Decontract P matrix block
     324        18568 :                   pmat = 0.0_dp
     325          460 :                   DO i = 1, SIZE(ec_env%matrix_p, 1)
     326          230 :                      IF (i == 1) THEN
     327          230 :                         pblock => p_block_a
     328              :                      ELSE
     329            0 :                         pblock => p_block_b
     330              :                      END IF
     331          230 :                      IF (.NOT. ASSOCIATED(pblock)) CYCLE
     332          230 :                      IF (trans) THEN
     333              :                         CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), &
     334              :                                    1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     335              :                                    pblock(sgfb, sgfa), SIZE(pblock, 1), &
     336           78 :                                    0.0_dp, work(1, 1), ldab)
     337              :                      ELSE
     338              :                         CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
     339              :                                    1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
     340              :                                    pblock(sgfa, sgfb), SIZE(pblock, 1), &
     341          152 :                                    0.0_dp, work(1, 1), ldab)
     342              :                      END IF
     343              :                      CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
     344              :                                 1.0_dp, work(1, 1), ldab, &
     345              :                                 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
     346          460 :                                 1.0_dp, pmat(1, 1), ncoa)
     347              :                   END DO
     348              : 
     349              :                   CALL dipole_force(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
     350              :                                     lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
     351          230 :                                     1, rac, rbc, pmat, forcea, forceb)
     352              : 
     353          462 :                   DEALLOCATE (work, pmat)
     354              :                END DO
     355              :             END DO
     356              : 
     357          322 :             DO idir = 1, 3
     358              :                force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) &
     359          924 :                                                   + fdir*fieldpol(idir)*forcea(idir, 1:3)
     360              :                force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) &
     361         1001 :                                                   + fdir*fieldpol(idir)*forceb(idir, 1:3)
     362              :             END DO
     363              : 
     364              :          END DO
     365           14 :          CALL neighbor_list_iterator_release(nl_iterator)
     366           14 :          DEALLOCATE (basis_set_list)
     367              :       END IF
     368              : 
     369           86 :       CALL timestop(handle)
     370              : 
     371          172 :    END SUBROUTINE ec_efield_mo_derivatives
     372              : 
     373              : END MODULE ec_efield_local
        

Generated by: LCOV version 2.0-1