LCOV - code coverage report
Current view: top level - src - qs_dcdr_ao.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 237 246 96.3 %
Date: 2024-03-28 07:31:50 Functions: 6 6 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 Calculate the derivatives of the MO coefficients wrt nuclear coordinates
      10             : !> \author Sandra Luber, Edward Ditler
      11             : ! **************************************************************************************************
      12             : 
      13             : MODULE qs_dcdr_ao
      14             : 
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      16             :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      17             :                                               gto_basis_set_type
      18             :    USE core_ae,                         ONLY: build_core_ae
      19             :    USE core_ppl,                        ONLY: build_core_ppl
      20             :    USE core_ppnl,                       ONLY: build_core_ppnl
      21             :    USE cp_control_types,                ONLY: dft_control_type
      22             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      23             :                                               copy_fm_to_dbcsr
      24             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25             :                                               cp_fm_release,&
      26             :                                               cp_fm_type
      27             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      28             :                                               cp_logger_type
      29             :    USE dbcsr_api,                       ONLY: dbcsr_copy,&
      30             :                                               dbcsr_get_block_p,&
      31             :                                               dbcsr_p_type,&
      32             :                                               dbcsr_set,&
      33             :                                               dbcsr_type
      34             :    USE input_constants,                 ONLY: do_ppl_analytic
      35             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      36             :                                               section_vals_type
      37             :    USE kinds,                           ONLY: default_string_length,&
      38             :                                               dp
      39             :    USE orbital_pointers,                ONLY: ncoset
      40             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      41             :    USE particle_types,                  ONLY: particle_type
      42             :    USE pw_env_types,                    ONLY: pw_env_get,&
      43             :                                               pw_env_type
      44             :    USE pw_methods,                      ONLY: pw_axpy,&
      45             :                                               pw_copy,&
      46             :                                               pw_scale,&
      47             :                                               pw_transfer,&
      48             :                                               pw_zero
      49             :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      50             :    USE pw_poisson_types,                ONLY: pw_poisson_type
      51             :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      52             :                                               pw_pool_type
      53             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      54             :                                               pw_r3d_rs_type
      55             :    USE qs_collocate_density,            ONLY: calculate_drho_core,&
      56             :                                               calculate_drho_elec_dR
      57             :    USE qs_energy_types,                 ONLY: qs_energy_type
      58             :    USE qs_environment_types,            ONLY: get_qs_env,&
      59             :                                               qs_environment_type
      60             :    USE qs_force_types,                  ONLY: qs_force_type
      61             :    USE qs_integral_utils,               ONLY: basis_set_list_setup,&
      62             :                                               get_memory_usage
      63             :    USE qs_integrate_potential,          ONLY: integrate_v_dbasis,&
      64             :                                               integrate_v_rspace
      65             :    USE qs_kind_types,                   ONLY: qs_kind_type
      66             :    USE qs_ks_types,                     ONLY: get_ks_env,&
      67             :                                               qs_ks_env_type
      68             :    USE qs_linres_types,                 ONLY: dcdr_env_type
      69             :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      70             :                                               get_neighbor_list_set_p,&
      71             :                                               neighbor_list_iterate,&
      72             :                                               neighbor_list_iterator_create,&
      73             :                                               neighbor_list_iterator_p_type,&
      74             :                                               neighbor_list_iterator_release,&
      75             :                                               neighbor_list_set_p_type
      76             :    USE qs_rho_methods,                  ONLY: qs_rho_rebuild,&
      77             :                                               qs_rho_update_rho
      78             :    USE qs_rho_types,                    ONLY: qs_rho_create,&
      79             :                                               qs_rho_get,&
      80             :                                               qs_rho_release,&
      81             :                                               qs_rho_type
      82             :    USE qs_vxc,                          ONLY: qs_vxc_create
      83             :    USE virial_types,                    ONLY: virial_type
      84             :    USE xc,                              ONLY: xc_calc_2nd_deriv,&
      85             :                                               xc_prep_2nd_deriv
      86             :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      87             :                                               xc_dset_release
      88             :    USE xc_rho_set_types,                ONLY: xc_rho_set_release,&
      89             :                                               xc_rho_set_type
      90             : 
      91             : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      92             : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
      93             : !$                    omp_init_lock, omp_set_lock, &
      94             : !$                    omp_unset_lock, omp_destroy_lock
      95             : 
      96             : #include "./base/base_uses.f90"
      97             : 
      98             :    IMPLICIT NONE
      99             : 
     100             :    PRIVATE
     101             :    PUBLIC :: core_dR, d_vhxc_dR, d_core_charge_density_dR, apply_op_constant_term
     102             :    PUBLIC :: vhxc_R_perturbed_basis_functions
     103             :    PUBLIC :: hr_mult_by_delta_1d
     104             : 
     105             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_ao'
     106             : 
     107             : CONTAINS
     108             : 
     109             : ! **************************************************************************************************
     110             : !> \brief Build the perturbed density matrix correction depending on the overlap derivative
     111             : !> \param qs_env ...
     112             : !> \param dcdr_env ...
     113             : !> \param overlap1 Overlap derivative in AO basis
     114             : !> \author Edward Ditler
     115             : ! **************************************************************************************************
     116         126 :    SUBROUTINE apply_op_constant_term(qs_env, dcdr_env, overlap1)
     117             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     118             :       TYPE(dcdr_env_type)                                :: dcdr_env
     119             :       TYPE(dbcsr_p_type), OPTIONAL                       :: overlap1
     120             : 
     121             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_constant_term'
     122             : 
     123             :       INTEGER                                            :: handle, ispin
     124             :       REAL(KIND=dp)                                      :: energy_hartree
     125             :       TYPE(cp_fm_type)                                   :: rho_ao_fm, rho_ao_s1, rho_ao_s1_rho_ao, &
     126             :                                                             s1_ao
     127         126 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho1_ao, rho_ao
     128             :       TYPE(pw_c1d_gs_type)                               :: rho1_tot_gspace, v_hartree_gspace
     129         126 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_g, rho1_g_pw
     130             :       TYPE(pw_env_type), POINTER                         :: pw_env
     131             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     132             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     133             :       TYPE(pw_r3d_rs_type)                               :: v_hartree_rspace
     134         252 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, rho_r, tau1_r, v_rspace_new, &
     135         126 :                                                             v_xc, v_xc_tau
     136             :       TYPE(qs_rho_type), POINTER                         :: perturbed_density, rho
     137             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     138             :       TYPE(xc_derivative_set_type)                       :: deriv_set
     139             :       TYPE(xc_rho_set_type)                              :: rho_set
     140             : 
     141             :       ! Build the perturbed density matrix correction depending on the overlap derivative
     142             :       !   P1 = C0 C1 + C1 C0
     143             :       !        - C0_(mu j) S1_(jk) C0_(k nu)
     144             :       ! This routine is adapted from apply_op_2_dft. There, build_dm_response builds
     145             :       !  C0 * dCR + dCR * C0.
     146             :       ! build_dm_response is computing $-1 * (C^0 C^1 + C^1 C^0)$ and later on in the
     147             :       !  integration the factor 2 is applied to account for the occupancy.
     148             :       ! The sign is negative because the kernel is on the RHS of the Sternheimer equation.
     149             :       !
     150             :       ! The correction factor in this routine needs to have
     151             :       !      the opposite sign mathematically as (C0 C1 + C1 C0)
     152             :       !   so the same sign in the code     because of the $-1$ in dCR
     153             :       !   so the opposite sign in the code because we are on the LHS of the Sternheimer equation.
     154             :       !
     155             :       ! This term must not go into the kernel applied by the linear response solver, because
     156             :       !  for the (P)CG algorithm, all constant terms have to be on one side of the equations
     157             :       !  and all solution dependent terms must be on the other side.
     158             : 
     159         126 :       CALL timeset(routineN, handle)
     160             : 
     161         126 :       NULLIFY (auxbas_pw_pool, pw_env, rho1_r, rho1_g_pw, &
     162         126 :                v_xc, poisson_env, input, rho, rho1_g, v_xc_tau)
     163             : 
     164         126 :       CALL cp_fm_create(rho_ao_fm, dcdr_env%aoao_fm_struct)
     165         126 :       CALL cp_fm_create(rho_ao_s1, dcdr_env%aoao_fm_struct)
     166         126 :       CALL cp_fm_create(rho_ao_s1_rho_ao, dcdr_env%aoao_fm_struct)
     167         126 :       CALL cp_fm_create(s1_ao, dcdr_env%aoao_fm_struct)
     168             : 
     169         126 :       IF (PRESENT(overlap1)) THEN
     170           0 :          CALL copy_dbcsr_to_fm(overlap1%matrix, s1_ao)
     171             :       ELSE
     172         126 :          CALL copy_dbcsr_to_fm(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, s1_ao)
     173             :       END IF
     174             : 
     175         288 :       DO ispin = 1, dcdr_env%nspins
     176         162 :          CALL dbcsr_set(dcdr_env%perturbed_dm_correction(ispin)%matrix, 0._dp)
     177         162 :          CALL dbcsr_set(dcdr_env%matrix_apply_op_constant(ispin)%matrix, 0.0_dp)
     178             : 
     179             :          CALL parallel_gemm('N', 'T', dcdr_env%nao, dcdr_env%nao, dcdr_env%nmo(ispin), &
     180             :                             1.0_dp, dcdr_env%mo_coeff(ispin), dcdr_env%mo_coeff(ispin), &
     181         162 :                             0.0_dp, rho_ao_fm)
     182             : 
     183             :          CALL parallel_gemm('N', 'N', dcdr_env%nao, dcdr_env%nao, dcdr_env%nao, &
     184             :                             1.0_dp, rho_ao_fm, s1_ao, &
     185         162 :                             0.0_dp, rho_ao_s1)
     186             : 
     187             :          CALL parallel_gemm('N', 'N', dcdr_env%nao, dcdr_env%nao, dcdr_env%nao, &
     188             :                             -1._dp, rho_ao_s1, rho_ao_fm, &   ! this is the sign mentioned above.
     189         162 :                             0.0_dp, rho_ao_s1_rho_ao)
     190             : 
     191         288 :          CALL copy_fm_to_dbcsr(rho_ao_s1_rho_ao, dcdr_env%perturbed_dm_correction(ispin)%matrix)
     192             :       END DO
     193             : 
     194         126 :       CALL cp_fm_release(rho_ao_fm)
     195         126 :       CALL cp_fm_release(rho_ao_s1)
     196         126 :       CALL cp_fm_release(rho_ao_s1_rho_ao)
     197         126 :       CALL cp_fm_release(s1_ao)
     198             :       ! Done building the density matrix correction
     199             : 
     200             :       ! Build the density struct from the environment
     201         126 :       NULLIFY (perturbed_density)
     202         126 :       ALLOCATE (perturbed_density)
     203         126 :       CALL qs_rho_create(perturbed_density)
     204         126 :       CALL qs_rho_rebuild(perturbed_density, qs_env=qs_env)
     205             : 
     206             :       ! ... set the density matrix to be the perturbed density matrix
     207         126 :       CALL qs_rho_get(perturbed_density, rho_ao=rho1_ao)
     208         288 :       DO ispin = 1, dcdr_env%nspins
     209         288 :          CALL dbcsr_copy(rho1_ao(ispin)%matrix, dcdr_env%perturbed_dm_correction(ispin)%matrix)
     210             :       END DO
     211             : 
     212             :       ! ... updates rho_r and rho_g to the rho%rho_ao.
     213             :       CALL qs_rho_update_rho(rho_struct=perturbed_density, &
     214         126 :                              qs_env=qs_env)
     215             : 
     216             :       ! Also update the qs_env%rho
     217         126 :       CALL get_qs_env(qs_env, rho=rho)
     218         126 :       CALL qs_rho_update_rho(rho, qs_env=qs_env)
     219         126 :       CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
     220             : 
     221         126 :       energy_hartree = 0.0_dp
     222             : 
     223             :       CALL get_qs_env(qs_env=qs_env, &
     224             :                       pw_env=pw_env, &
     225         126 :                       input=input)
     226             : 
     227             :       ! Create the temporary grids
     228             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     229         126 :                       poisson_env=poisson_env)
     230             : 
     231             :       ! Allocate deriv_set and rho_set
     232         126 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     233             : 
     234             :       CALL xc_prep_2nd_deriv(deriv_set, rho_set, &
     235             :                              rho_r, auxbas_pw_pool, &
     236         126 :                              xc_section=xc_section)
     237             : 
     238             :       ! Done with deriv_set and rho_set
     239             : 
     240         540 :       ALLOCATE (v_rspace_new(dcdr_env%nspins))
     241         126 :       CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
     242         126 :       CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
     243             : 
     244             :       ! Calculate the Hartree potential on the total density
     245         126 :       CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
     246             : 
     247         126 :       CALL qs_rho_get(perturbed_density, rho_g=rho1_g, rho_r=rho1_r, tau_r=tau1_r)
     248         126 :       CALL pw_copy(rho1_g(1), rho1_tot_gspace)
     249         162 :       DO ispin = 2, dcdr_env%nspins
     250         162 :          CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
     251             :       END DO
     252             : 
     253             :       CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
     254             :                             energy_hartree, &
     255         126 :                             v_hartree_gspace)
     256         126 :       CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
     257             : 
     258         126 :       CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
     259             : 
     260             :       ! Calculate the second derivative of the exchange-correlation potential
     261             :       CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, &
     262         126 :                              rho1_r, rho1_g_pw, tau1_r, auxbas_pw_pool, xc_section, gapw=.FALSE.)
     263             : 
     264         288 :       DO ispin = 1, dcdr_env%nspins
     265         288 :          v_rspace_new(ispin) = v_xc(ispin)
     266             :       END DO
     267         126 :       DEALLOCATE (v_xc)
     268             : 
     269             :       ! Done calculating the potentials
     270             : 
     271             :       !-------------------------------!
     272             :       ! Add both hartree and xc terms !
     273             :       !-------------------------------!
     274         126 :       CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
     275         288 :       DO ispin = 1, dcdr_env%nspins
     276         288 :          CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
     277             :       END DO
     278             : 
     279         288 :       DO ispin = 1, dcdr_env%nspins
     280         162 :          CALL dbcsr_set(dcdr_env%matrix_apply_op_constant(ispin)%matrix, 0.0_dp)
     281         162 :          CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin))
     282         162 :          IF (dcdr_env%nspins == 1) THEN
     283          90 :             CALL pw_scale(v_rspace_new(1), 2.0_dp)
     284             :          END IF
     285             : 
     286             :          CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
     287             :                                  hmat=dcdr_env%matrix_apply_op_constant(ispin), &
     288             :                                  qs_env=qs_env, &
     289         288 :                                  calculate_forces=.FALSE.)
     290             :       END DO
     291             : 
     292         126 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
     293         126 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
     294         288 :       DO ispin = 1, dcdr_env%nspins
     295         288 :          CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
     296             :       END DO
     297         126 :       DEALLOCATE (v_rspace_new)
     298             : 
     299         126 :       IF (ASSOCIATED(v_xc_tau)) THEN
     300           0 :          CALL pw_scale(v_xc_tau(1), 2._dp*v_xc_tau(1)%pw_grid%dvol)
     301             :          CALL integrate_v_rspace(v_rspace=v_xc_tau(1), &
     302             :                                  hmat=dcdr_env%matrix_apply_op_constant(1), &
     303             :                                  qs_env=qs_env, &
     304             :                                  compute_tau=.TRUE., &
     305           0 :                                  calculate_forces=.FALSE.)
     306             : 
     307           0 :          CALL auxbas_pw_pool%give_back_pw(v_xc_tau(1))
     308           0 :          DEALLOCATE (v_xc_tau)
     309             :       END IF
     310             : 
     311         126 :       CALL qs_rho_release(perturbed_density)
     312         126 :       DEALLOCATE (perturbed_density)
     313         126 :       CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
     314         126 :       CALL xc_dset_release(deriv_set)
     315             : 
     316         126 :       CALL timestop(handle)
     317             : 
     318        2898 :    END SUBROUTINE apply_op_constant_term
     319             : 
     320             : ! **************************************************************************************************
     321             : !> \brief Calculate the derivative of the Hartree term due to the core charge density
     322             : !> \param qs_env ...
     323             : !> \param dcdr_env ...
     324             : !> \author Edward Ditler
     325             : ! **************************************************************************************************
     326          36 :    SUBROUTINE d_core_charge_density_dR(qs_env, dcdr_env)
     327             :       ! drho_core contribution
     328             :       ! sum over all directions
     329             :       ! output in ao x ao
     330             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     331             :       TYPE(dcdr_env_type)                                :: dcdr_env
     332             : 
     333             :       CHARACTER(len=*), PARAMETER :: routineN = 'd_core_charge_density_dR'
     334             : 
     335             :       INTEGER                                            :: beta, handle
     336             :       TYPE(cp_logger_type), POINTER                      :: logger
     337             :       TYPE(dft_control_type), POINTER                    :: dft_control
     338             :       TYPE(pw_c1d_gs_type)                               :: drho_g, v_hartree_gspace
     339             :       TYPE(pw_env_type), POINTER                         :: pw_env
     340             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     341          36 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     342             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     343             :       TYPE(pw_r3d_rs_type)                               :: v_hartree_rspace
     344             :       TYPE(qs_rho_type), POINTER                         :: rho
     345             : 
     346          36 :       CALL timeset(routineN, handle)
     347             : 
     348          36 :       logger => cp_get_default_logger()
     349             : 
     350          36 :       NULLIFY (pw_env, auxbas_pw_pool, pw_pools, poisson_env, dft_control, &
     351          36 :                rho)
     352             : 
     353             :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, &
     354          36 :                       dft_control=dft_control)
     355             : 
     356             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env, &
     357          36 :                       pw_pools=pw_pools)
     358             : 
     359             :       ! Create the Hartree potential grids in real and reciprocal space.
     360          36 :       CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
     361          36 :       CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
     362             :       ! Create the grid for the derivative of the core potential
     363          36 :       CALL auxbas_pw_pool%create_pw(drho_g)
     364             : 
     365         144 :       DO beta = 1, 3
     366         108 :          CALL pw_zero(v_hartree_gspace)
     367         108 :          CALL pw_zero(v_hartree_rspace)
     368         108 :          CALL pw_zero(drho_g)
     369             : 
     370             :          ! Calculate the Hartree potential on the perturbed density and Poisson solve it
     371             :          CALL calculate_drho_core(drho_core=drho_g, qs_env=qs_env, &
     372         108 :                                   beta=beta, lambda=dcdr_env%lambda)
     373             :          CALL pw_poisson_solve(poisson_env, drho_g, &
     374         108 :                                vhartree=v_hartree_gspace)
     375         108 :          CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
     376         108 :          CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
     377             : 
     378             :          ! Calculate the integrals
     379             :          CALL integrate_v_rspace(v_rspace=v_hartree_rspace, &
     380             :                                  hmat=dcdr_env%matrix_core_charge_1(beta), &
     381             :                                  qs_env=qs_env, &
     382         144 :                                  calculate_forces=.FALSE.)
     383             :       END DO
     384             : 
     385          36 :       CALL auxbas_pw_pool%give_back_pw(drho_g)
     386          36 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
     387          36 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
     388             : 
     389          36 :       CALL timestop(handle)
     390          36 :    END SUBROUTINE d_core_charge_density_dR
     391             : 
     392             : ! **************************************************************************************************
     393             : !> \brief Core Hamiltonian contributions to the operator (the pseudopotentials)
     394             : !> \param qs_env ...
     395             : !> \param dcdr_env ..
     396             : !> \author Edward Ditler
     397             : ! **************************************************************************************************
     398          36 :    SUBROUTINE core_dR(qs_env, dcdr_env)
     399             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     400             :       TYPE(dcdr_env_type)                                :: dcdr_env
     401             : 
     402             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'core_dR'
     403             : 
     404             :       CHARACTER(LEN=default_string_length)               :: my_basis_type
     405             :       INTEGER                                            :: handle, nder
     406          36 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     407             :       LOGICAL                                            :: calculate_forces, failure, ppl_present, &
     408             :                                                             ppnl_present, use_virial
     409             :       REAL(KIND=dp)                                      :: eps_ppnl
     410             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: deltaR
     411          36 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     412          36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
     413          36 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_hc_pass, matrix_p_pass, &
     414          36 :                                                             matrix_ppnl_1_pass
     415             :       TYPE(dft_control_type), POINTER                    :: dft_control
     416             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     417          36 :          POINTER                                         :: sab_orb, sac_ae, sac_ppl, sap_ppnl
     418          36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     419          36 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     420          36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     421             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     422             :       TYPE(qs_rho_type), POINTER                         :: rho
     423             :       TYPE(virial_type), POINTER                         :: virial
     424             : 
     425          36 :       CALL timeset(routineN, handle)
     426             : 
     427          36 :       failure = .FALSE.
     428             : 
     429          36 :       NULLIFY (atomic_kind_set, qs_kind_set, ks_env, dft_control, particle_set, &
     430          36 :                sab_orb, sac_ae, sac_ppl, sap_ppnl, virial, rho, rho_ao)
     431             : 
     432             :       CALL get_qs_env(qs_env=qs_env, &
     433             :                       atomic_kind_set=atomic_kind_set, &
     434             :                       qs_kind_set=qs_kind_set, &
     435             :                       ks_env=ks_env, &
     436             :                       dft_control=dft_control, &
     437             :                       particle_set=particle_set, &
     438             :                       sab_orb=sab_orb, &
     439             :                       sac_ae=sac_ae, &
     440             :                       sac_ppl=sac_ppl, &
     441             :                       sap_ppnl=sap_ppnl, &
     442          36 :                       virial=virial)
     443          36 :       CALL get_ks_env(ks_env=ks_env, rho=rho)
     444          36 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     445          36 :       deltaR => dcdr_env%delta_basis_function
     446             : 
     447          36 :       nder = 1
     448          36 :       calculate_forces = .FALSE.
     449             : 
     450          36 :       my_basis_type = "ORB"
     451             : 
     452             :       ! ECP/AE contribution to the core hamiltonian
     453          36 :       IF (ASSOCIATED(sac_ae)) THEN
     454           0 :          CPABORT("ECP/AE functionality in qs_dcdr_ao missing")
     455             :          ! Missing feature: deltaR weighting factors of the derivatives wrt. nuclear positions
     456           0 :          matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
     457           0 :          matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
     458             :          CALL build_core_ae(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, &
     459             :                             force=force, virial=virial, calculate_forces=calculate_forces, &
     460             :                             use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, &
     461             :                             atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
     462             :                             sab_orb=sab_orb, sac_ae=sac_ae, nimages=1, &
     463           0 :                             cell_to_index=cell_to_index)
     464             :       END IF
     465             :       ! *** compute the ppl contribution to the core hamiltonian ***
     466          36 :       ppl_present = ASSOCIATED(sac_ppl)
     467          36 :       IF (ppl_present) THEN
     468          36 :          IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
     469          36 :             matrix_hc_pass(1:3, 1:1) => dcdr_env%matrix_hc(1:3)
     470          36 :             matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
     471             :             CALL build_core_ppl(matrix_h=matrix_hc_pass, matrix_p=matrix_p_pass, &
     472             :                                 force=force, virial=virial, calculate_forces=calculate_forces, &
     473             :                                 use_virial=use_virial, nder=nder, qs_kind_set=qs_kind_set, &
     474             :                                 atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
     475             :                                 sab_orb=sab_orb, sac_ppl=sac_ppl, basis_type=my_basis_type, &
     476          36 :                                 nimages=1, cell_to_index=cell_to_index, deltaR=deltaR)
     477             : 
     478             :          END IF ! ppl_analytic
     479             :       END IF ! ppl_present
     480             : 
     481             :       ! *** compute the ppnl contribution to the core hamiltonian ***
     482          36 :       eps_ppnl = dft_control%qs_control%eps_ppnl
     483          36 :       ppnl_present = ASSOCIATED(sap_ppnl)
     484          36 :       IF (ppnl_present) THEN
     485          36 :          matrix_ppnl_1_pass(1:3, 1:1) => dcdr_env%matrix_ppnl_1(1:3)
     486             :          CALL build_core_ppnl(matrix_h=matrix_ppnl_1_pass, matrix_p=matrix_p_pass, force=force, virial=virial, &
     487             :                               calculate_forces=calculate_forces, use_virial=use_virial, nder=nder, &
     488             :                               qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
     489             :                               particle_set=particle_set, sab_orb=sab_orb, sap_ppnl=sap_ppnl, &
     490             :                               eps_ppnl=eps_ppnl, nimages=1, cell_to_index=cell_to_index, &
     491          36 :                               basis_type=my_basis_type, deltaR=deltaR)
     492             :       END IF
     493             : 
     494          36 :       CALL timestop(handle)
     495          36 :    END SUBROUTINE core_dR
     496             : 
     497             : ! **************************************************************************************************
     498             : !> \brief The derivatives of the basis functions going into the HXC potential wrt nuclear positions
     499             : !> \param qs_env ...
     500             : !> \param dcdr_env ...
     501             : !> \author Edward Ditler
     502             : ! **************************************************************************************************
     503          36 :    SUBROUTINE d_vhxc_dR(qs_env, dcdr_env)
     504             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     505             :       TYPE(dcdr_env_type)                                :: dcdr_env
     506             : 
     507             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'd_vhxc_dR'
     508             : 
     509             :       INTEGER                                            :: handle, idir, ispin
     510          36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
     511             :       TYPE(pw_c1d_gs_type)                               :: drho_g_total, v_hartree_gspace
     512          36 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: drho_g
     513             :       TYPE(pw_env_type), POINTER                         :: pw_env
     514             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     515          36 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     516             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     517             :       TYPE(pw_r3d_rs_type)                               :: drho_r_total, v_hartree_rspace
     518          72 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: drho_r, dtau_r, rho_r, v_xc, v_xc_tau
     519             :       TYPE(qs_rho_type), POINTER                         :: rho
     520             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     521             :       TYPE(xc_derivative_set_type)                       :: my_deriv_set
     522             :       TYPE(xc_rho_set_type)                              :: my_rho_set
     523             : 
     524          36 :       CALL timeset(routineN, handle)
     525             : 
     526             :       CALL get_qs_env(qs_env=qs_env, &
     527             :                       pw_env=pw_env, &
     528             :                       input=input, &
     529          36 :                       rho=rho)
     530          36 :       CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
     531             : 
     532          36 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     533             : 
     534             :       ! get the tmp grids
     535         150 :       ALLOCATE (drho_r(dcdr_env%nspins))
     536         150 :       ALLOCATE (drho_g(dcdr_env%nspins))
     537             : 
     538             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     539          36 :                       pw_pools=pw_pools, poisson_env=poisson_env)
     540          36 :       CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
     541          36 :       CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
     542             : 
     543          78 :       DO ispin = 1, dcdr_env%nspins
     544          42 :          CALL auxbas_pw_pool%create_pw(drho_r(ispin))
     545          78 :          CALL auxbas_pw_pool%create_pw(drho_g(ispin))
     546             :       END DO
     547          36 :       CALL auxbas_pw_pool%create_pw(drho_g_total)
     548          36 :       CALL auxbas_pw_pool%create_pw(drho_r_total)
     549             : 
     550         144 :       DO idir = 1, 3
     551         108 :          CALL pw_zero(v_hartree_gspace)
     552         108 :          CALL pw_zero(v_hartree_rspace)
     553         108 :          CALL pw_zero(drho_g_total)
     554         108 :          CALL pw_zero(drho_r_total)
     555             : 
     556         234 :          DO ispin = 1, dcdr_env%nspins
     557         126 :             CALL pw_zero(drho_r(ispin))
     558         126 :             CALL pw_zero(drho_g(ispin))
     559             : 
     560             :             ! Get the density
     561             :             CALL calculate_drho_elec_dR(matrix_p=rho_ao(ispin)%matrix, &
     562             :                                         drho=drho_r(ispin), &
     563             :                                         drho_gspace=drho_g(ispin), &
     564             :                                         qs_env=qs_env, &
     565         126 :                                         beta=idir, lambda=dcdr_env%lambda)
     566             : 
     567         126 :             CALL pw_axpy(drho_g(ispin), drho_g_total)
     568         234 :             CALL pw_axpy(drho_r(ispin), drho_r_total)
     569             :          END DO
     570             :          ! Get the Hartree potential corresponding to the perturbed density
     571             :          CALL pw_poisson_solve(poisson_env, drho_g_total, &
     572         108 :                                vhartree=v_hartree_gspace)
     573         108 :          CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
     574             : 
     575             :          ! Get the XC potential corresponding to the perturbed density
     576             :          CALL xc_prep_2nd_deriv(my_deriv_set, my_rho_set, &
     577             :                                 rho_r, auxbas_pw_pool, &
     578         108 :                                 xc_section=xc_section)
     579             : 
     580         108 :          NULLIFY (dtau_r)
     581             :          CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, my_deriv_set, my_rho_set, &
     582         108 :                                 drho_r, drho_g, dtau_r, auxbas_pw_pool, xc_section, gapw=.FALSE.)
     583         108 :          IF (ASSOCIATED(v_xc_tau)) CPABORT("Meta functionals are not supported!")
     584             : 
     585         108 :          CALL xc_dset_release(my_deriv_set)
     586         108 :          CALL xc_rho_set_release(my_rho_set)
     587             : 
     588             :          !-------------------------------!
     589             :          ! Add both hartree and xc terms !
     590             :          !-------------------------------!
     591         234 :          DO ispin = 1, dcdr_env%nspins
     592             :             ! Can the dvol be different?
     593         126 :             CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
     594         126 :             CALL pw_axpy(v_hartree_rspace, v_xc(ispin), v_hartree_rspace%pw_grid%dvol)
     595             : 
     596             :             CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
     597             :                                     hmat=dcdr_env%matrix_d_vhxc_dR(idir, ispin), &
     598             :                                     qs_env=qs_env, &
     599         126 :                                     calculate_forces=.FALSE.)
     600             : 
     601             :             ! v_xc gets allocated again in xc_calc_2nd_deriv
     602         234 :             CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
     603             :          END DO ! ispin
     604         144 :          DEALLOCATE (v_xc)
     605             :       END DO ! idir
     606             : 
     607          36 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
     608          36 :       CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
     609          36 :       CALL auxbas_pw_pool%give_back_pw(drho_g_total)
     610          36 :       CALL auxbas_pw_pool%give_back_pw(drho_r_total)
     611             : 
     612          78 :       DO ispin = 1, dcdr_env%nspins
     613          42 :          CALL auxbas_pw_pool%give_back_pw(drho_g(ispin))
     614          78 :          CALL auxbas_pw_pool%give_back_pw(drho_r(ispin))
     615             :       END DO
     616             : 
     617          36 :       DEALLOCATE (drho_g)
     618          36 :       DEALLOCATE (drho_r)
     619             : 
     620          36 :       CALL timestop(handle)
     621             : 
     622         792 :    END SUBROUTINE d_vhxc_dR
     623             : 
     624             : ! **************************************************************************************************
     625             : !> \brief The derivatives of the basis functions over which the HXC potential is integrated,
     626             : !>          so < da/dR | Vhxc | b >
     627             : !> \param qs_env ...
     628             : !> \param dcdr_env ...
     629             : !> \author Edward Ditler
     630             : ! **************************************************************************************************
     631          36 :    SUBROUTINE vhxc_R_perturbed_basis_functions(qs_env, dcdr_env)
     632             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     633             :       TYPE(dcdr_env_type)                                :: dcdr_env
     634             : 
     635             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'vhxc_R_perturbed_basis_functions'
     636             : 
     637             :       INTEGER                                            :: handle, ispin
     638          36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_vhxc_dbasis
     639          36 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p
     640             :       TYPE(pw_env_type), POINTER                         :: pw_env
     641             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     642          36 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_hxc_r, v_tau_rspace
     643             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_r
     644             :       TYPE(qs_energy_type), POINTER                      :: energy
     645             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     646             :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     647             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     648             : 
     649          36 :       CALL timeset(routineN, handle)
     650             : 
     651          36 :       NULLIFY (rho_struct, energy, input, ks_env, pw_env, matrix_p)
     652             :       CALL get_qs_env(qs_env, &
     653             :                       rho=rho_struct, &
     654             :                       energy=energy, &
     655             :                       input=input, &
     656             :                       ks_env=ks_env, &
     657             :                       pw_env=pw_env, &
     658          36 :                       v_hartree_rspace=v_hartree_r)
     659          36 :       CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
     660          36 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     661             : 
     662          36 :       NULLIFY (auxbas_pw_pool)
     663          36 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     664             : 
     665             :       ! *** calculate the xc potential on the pw density ***
     666             :       ! *** associates v_hxc_r if the xc potential needs to be computed.
     667             :       ! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set
     668          36 :       NULLIFY (v_hxc_r, v_tau_rspace)
     669             :       CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
     670          36 :                          vxc_rho=v_hxc_r, vxc_tau=v_tau_rspace, exc=energy%exc)
     671             : 
     672          78 :       DO ispin = 1, dcdr_env%nspins
     673          42 :          CALL pw_scale(v_hxc_r(ispin), v_hxc_r(ispin)%pw_grid%dvol)
     674             : 
     675             :          ! sum up potentials and integrate
     676          42 :          CALL pw_axpy(v_hartree_r, v_hxc_r(ispin), 1._dp)
     677             : 
     678          42 :          matrix_vhxc_dbasis => dcdr_env%matrix_vhxc_perturbed_basis(ispin, :)
     679             :          CALL integrate_v_dbasis(v_rspace=v_hxc_r(ispin), &
     680             :                                  matrix_p=matrix_p(ispin, 1)%matrix, &
     681             :                                  matrix_vhxc_dbasis=matrix_vhxc_dbasis, &
     682             :                                  qs_env=qs_env, &
     683          42 :                                  lambda=dcdr_env%lambda)
     684             : 
     685          78 :          CALL auxbas_pw_pool%give_back_pw(v_hxc_r(ispin))
     686             :       END DO
     687             : 
     688          36 :       DEALLOCATE (v_hxc_r)
     689             : 
     690          36 :       CALL timestop(handle)
     691          36 :    END SUBROUTINE vhxc_R_perturbed_basis_functions
     692             : 
     693             : ! **************************************************************************************************
     694             : !> \brief Enforce that one of the basis functions in < a | O | b > is centered on atom lambda.
     695             : !> \param matrix ...
     696             : !> \param qs_kind_set ...
     697             : !> \param basis_type ...
     698             : !> \param sab_nl ...
     699             : !> \param lambda Atom index
     700             : !> \param direction_Or True: < a | O | b==lambda >, False: < a==lambda | O | b >
     701             : ! **************************************************************************************************
     702        2178 :    SUBROUTINE hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_Or)
     703             : 
     704             :       TYPE(dbcsr_type), POINTER                          :: matrix
     705             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     706             :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
     707             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     708             :          POINTER                                         :: sab_nl
     709             :       INTEGER, INTENT(IN)                                :: lambda
     710             :       LOGICAL, INTENT(IN)                                :: direction_Or
     711             : 
     712             :       CHARACTER(len=*), PARAMETER :: routineN = 'hr_mult_by_delta_1d'
     713             : 
     714             :       INTEGER                                            :: handle, iatom, icol, ikind, irow, jatom, &
     715             :                                                             jkind, ldsab, mepos, nkind, nseta, &
     716             :                                                             nsetb, nthread
     717             :       INTEGER, DIMENSION(3)                              :: cell
     718        2178 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     719        2178 :                                                             npgfb, nsgfa, nsgfb
     720        2178 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     721             :       LOGICAL                                            :: do_symmetric, found
     722             :       REAL(KIND=dp), DIMENSION(3)                        :: rab
     723        2178 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     724        2178 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: k_block, rpgfa, rpgfb, scon_a, scon_b, &
     725        2178 :                                                             zeta, zetb
     726        2178 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     727             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     728             :       TYPE(neighbor_list_iterator_p_type), &
     729        2178 :          DIMENSION(:), POINTER                           :: nl_iterator
     730             : 
     731        2178 :       CALL timeset(routineN, handle)
     732             : 
     733        2178 :       nkind = SIZE(qs_kind_set)
     734             : 
     735             :       ! check for symmetry
     736        2178 :       CPASSERT(SIZE(sab_nl) > 0)
     737        2178 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
     738             : 
     739             :       ! prepare basis set
     740       10890 :       ALLOCATE (basis_set_list(nkind))
     741        2178 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
     742             : 
     743             :       ! *** Allocate work storage ***
     744        2178 :       ldsab = get_memory_usage(qs_kind_set, basis_type)
     745             : 
     746             :       nthread = 1
     747        2178 : !$    nthread = omp_get_max_threads()
     748             :       ! Iterate of neighbor list
     749        2178 :       CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
     750             : 
     751             : !$OMP PARALLEL DEFAULT(NONE) &
     752             : !$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric) &
     753             : !$OMP SHARED (ncoset,matrix,basis_set_list) &
     754             : !$OMP SHARED (direction_or, lambda) &
     755             : !$OMP PRIVATE (k_block,mepos,ikind,jkind,iatom,jatom,rab,cell) &
     756             : !$OMP PRIVATE (basis_set_a,basis_set_b) &
     757             : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a) &
     758             : !$OMP PRIVATE (zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb) &
     759        2178 : !$OMP PRIVATE (zetb, scon_a, scon_b, irow, icol, found)
     760             : 
     761             :       mepos = 0
     762             : !$    mepos = omp_get_thread_num()
     763             : 
     764             :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     765             :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
     766             :                                 iatom=iatom, jatom=jatom, r=rab, cell=cell)
     767             :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     768             :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     769             :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     770             :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     771             :          ! basis ikind
     772             :          first_sgfa => basis_set_a%first_sgf
     773             :          la_max => basis_set_a%lmax
     774             :          la_min => basis_set_a%lmin
     775             :          npgfa => basis_set_a%npgf
     776             :          nseta = basis_set_a%nset
     777             :          nsgfa => basis_set_a%nsgf_set
     778             :          rpgfa => basis_set_a%pgf_radius
     779             :          set_radius_a => basis_set_a%set_radius
     780             :          scon_a => basis_set_a%scon
     781             :          zeta => basis_set_a%zet
     782             :          ! basis jkind
     783             :          first_sgfb => basis_set_b%first_sgf
     784             :          lb_max => basis_set_b%lmax
     785             :          lb_min => basis_set_b%lmin
     786             :          npgfb => basis_set_b%npgf
     787             :          nsetb = basis_set_b%nset
     788             :          nsgfb => basis_set_b%nsgf_set
     789             :          rpgfb => basis_set_b%pgf_radius
     790             :          set_radius_b => basis_set_b%set_radius
     791             :          scon_b => basis_set_b%scon
     792             :          zetb => basis_set_b%zet
     793             : 
     794             :          IF (do_symmetric) THEN
     795             :             IF (iatom <= jatom) THEN
     796             :                irow = iatom
     797             :                icol = jatom
     798             :             ELSE
     799             :                irow = jatom
     800             :                icol = iatom
     801             :             END IF
     802             :          ELSE
     803             :             irow = iatom
     804             :             icol = jatom
     805             :          END IF
     806             : 
     807             :          NULLIFY (k_block)
     808             :          CALL dbcsr_get_block_p(matrix, irow, icol, k_block, found)
     809             :          CPASSERT(found)
     810             : 
     811             :          IF (direction_Or) THEN
     812             :             IF (jatom /= lambda) k_block(:, :) = 0._dp
     813             :          ELSE IF (.NOT. direction_Or) THEN
     814             :             IF (iatom /= lambda) k_block(:, :) = 0._dp
     815             :          END IF
     816             :       END DO
     817             : !$OMP END PARALLEL
     818        2178 :       CALL neighbor_list_iterator_release(nl_iterator)
     819             : 
     820             :       ! Release work storage
     821        2178 :       DEALLOCATE (basis_set_list)
     822             : 
     823        2178 :       CALL timestop(handle)
     824             : 
     825        4356 :    END SUBROUTINE hr_mult_by_delta_1d
     826             : 
     827             : END MODULE qs_dcdr_ao

Generated by: LCOV version 1.15