LCOV - code coverage report
Current view: top level - src - qs_dcdr.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.3 % 299 291
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 5 5

            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 Calculate the derivatives of the MO coefficients wrt nuclear coordinates
      10              : !> \author Sandra Luber, Edward Ditler
      11              : ! **************************************************************************************************
      12              : 
      13              : MODULE qs_dcdr
      14              : 
      15              :    USE atomic_kind_types,               ONLY: get_atomic_kind
      16              :    USE cell_types,                      ONLY: cell_type,&
      17              :                                               pbc
      18              :    USE cp_array_utils,                  ONLY: cp_2d_r_p_type
      19              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      20              :                                               dbcsr_copy,&
      21              :                                               dbcsr_desymmetrize,&
      22              :                                               dbcsr_p_type,&
      23              :                                               dbcsr_set
      24              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      25              :                                               dbcsr_allocate_matrix_set,&
      26              :                                               dbcsr_deallocate_matrix_set
      27              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      28              :                                               cp_fm_scale_and_add,&
      29              :                                               cp_fm_trace
      30              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      31              :                                               cp_fm_get_diag,&
      32              :                                               cp_fm_release,&
      33              :                                               cp_fm_set_all,&
      34              :                                               cp_fm_to_fm,&
      35              :                                               cp_fm_type
      36              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      37              :                                               cp_logger_type
      38              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      39              :                                               cp_print_key_unit_nr
      40              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      41              :                                               section_vals_type
      42              :    USE kinds,                           ONLY: dp
      43              :    USE molecule_types,                  ONLY: molecule_of_atom,&
      44              :                                               molecule_type
      45              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      46              :    USE particle_types,                  ONLY: particle_type
      47              :    USE qs_dcdr_ao,                      ONLY: apply_op_constant_term,&
      48              :                                               core_dR,&
      49              :                                               d_core_charge_density_dR,&
      50              :                                               d_vhxc_dR,&
      51              :                                               hr_mult_by_delta_1d,&
      52              :                                               vhxc_R_perturbed_basis_functions
      53              :    USE qs_dcdr_utils,                   ONLY: dcdr_read_restart,&
      54              :                                               dcdr_write_restart,&
      55              :                                               multiply_localization,&
      56              :                                               shift_wannier_into_cell
      57              :    USE qs_environment_types,            ONLY: get_qs_env,&
      58              :                                               qs_environment_type
      59              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      60              :                                               qs_kind_type
      61              :    USE qs_linres_methods,               ONLY: linres_solver
      62              :    USE qs_linres_types,                 ONLY: dcdr_env_type,&
      63              :                                               get_polar_env,&
      64              :                                               linres_control_type,&
      65              :                                               polar_env_type
      66              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      67              :                                               mo_set_type
      68              :    USE qs_moments,                      ONLY: build_local_moment_matrix,&
      69              :                                               dipole_deriv_ao
      70              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      71              :    USE qs_p_env_types,                  ONLY: qs_p_env_type
      72              : #include "./base/base_uses.f90"
      73              : 
      74              :    IMPLICIT NONE
      75              : 
      76              :    PRIVATE
      77              :    PUBLIC :: prepare_per_atom, dcdr_response_dR, dcdr_build_op_dR, apt_dR, apt_dR_localization
      78              : 
      79              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr'
      80              : 
      81              : CONTAINS
      82              : 
      83              : ! **************************************************************************************************
      84              : !> \brief Prepare the environment for a choice of lambda
      85              : !> \param dcdr_env ...
      86              : !> \param qs_env ...
      87              : !> \author Edward Ditler
      88              : ! **************************************************************************************************
      89           72 :    SUBROUTINE prepare_per_atom(dcdr_env, qs_env)
      90              :       TYPE(dcdr_env_type)                                :: dcdr_env
      91              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      92              : 
      93              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'prepare_per_atom'
      94              : 
      95              :       INTEGER                                            :: handle, i, ispin, j, natom
      96              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      97           72 :          POINTER                                         :: sab_all
      98           72 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      99           72 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     100              : 
     101           72 :       CALL timeset(routineN, handle)
     102              : 
     103           72 :       NULLIFY (sab_all, qs_kind_set, particle_set)
     104              :       CALL get_qs_env(qs_env=qs_env, &
     105              :                       sab_all=sab_all, &
     106              :                       qs_kind_set=qs_kind_set, &
     107           72 :                       particle_set=particle_set)
     108              : 
     109           72 :       natom = SIZE(particle_set)
     110           72 :       IF (dcdr_env%distributed_origin) dcdr_env%ref_point(:) = particle_set(dcdr_env%lambda)%r(:)
     111              : 
     112          936 :       dcdr_env%delta_basis_function = 0._dp
     113          288 :       dcdr_env%delta_basis_function(:, dcdr_env%lambda) = 1._dp
     114              : 
     115              :       ! S matrix
     116              :       ! S1 = - < da/dr | b > * delta_a - < a | db/dr > * delta_b
     117              : 
     118              :       ! matrix_s(2:4) are anti-symmetric matrices and contain derivatives wrt. to < a |
     119              :       !               = < da/dR | b > = - < da/dr | b > = < a | db/dr >
     120              :       ! matrix_s1(2:4) = d/dR < a | b >
     121              :       !                and it's built as
     122              :       !                = - matrix_s      * delta_b  +  matrix_s      * delta_a
     123              :       !                = - < da/dR | b > * delta_b  +  < da/dR | b > * delta_a
     124              :       !                = + < da/dr | b > * delta_b  -  < da/dr | b > * delta_a
     125              :       !                = - < a | db/dr > * delta_b  -  < da/dr | b > * delta_a
     126              : 
     127          288 :       DO i = 1, 3
     128              :          ! S matrix
     129          216 :          CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
     130          216 :          CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_s1(1 + i)%matrix)
     131          216 :          CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
     132              : 
     133              :          CALL hr_mult_by_delta_1d(dcdr_env%matrix_s1(1 + i)%matrix, qs_kind_set, "ORB", &
     134          216 :                                   sab_all, dcdr_env%lambda, direction_Or=.TRUE.)
     135              :          CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
     136          216 :                                   sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
     137              : 
     138          216 :          CALL dbcsr_add(dcdr_env%matrix_s1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
     139          216 :          CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
     140              : 
     141              :          ! T matrix
     142          216 :          CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
     143          216 :          CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_t1(1 + i)%matrix)
     144          216 :          CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
     145              : 
     146              :          CALL hr_mult_by_delta_1d(dcdr_env%matrix_t1(1 + i)%matrix, qs_kind_set, "ORB", &
     147          216 :                                   sab_all, dcdr_env%lambda, direction_Or=.TRUE.)
     148              :          CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
     149          216 :                                   sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
     150              : 
     151          216 :          CALL dbcsr_add(dcdr_env%matrix_t1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
     152          288 :          CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
     153              :       END DO
     154              : 
     155              :       ! Operator:
     156          156 :       DO ispin = 1, dcdr_env%nspins
     157          408 :          DO i = 1, 3
     158          252 :             CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
     159          252 :             CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
     160          252 :             CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i)%matrix, 0.0_dp)
     161          252 :             CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i + 3)%matrix, 0.0_dp)
     162          252 :             CALL dbcsr_set(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, 0.0_dp)
     163          336 :             CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
     164              :          END DO
     165              :       END DO
     166              : 
     167           72 :       CALL core_dR(qs_env, dcdr_env)  ! dcdr_env%matrix_ppnl_1, hc
     168           72 :       CALL d_vhxc_dR(qs_env, dcdr_env)  ! dcdr_env%matrix_d_vhxc_dR
     169           72 :       CALL d_core_charge_density_dR(qs_env, dcdr_env)  ! dcdr_env%matrix_core_charge_1
     170           72 :       CALL vhxc_R_perturbed_basis_functions(qs_env, dcdr_env)  ! dcdr_env%matrix_vhxc_perturbed_basis
     171              : 
     172              :       ! APT:
     173          288 :       DO i = 1, 3
     174          936 :          DO j = 1, 3
     175          864 :             CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0._dp)
     176              :          END DO
     177              :       END DO
     178              : 
     179           72 :       CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, 1, dcdr_env%ref_point)
     180              : 
     181           72 :       CALL timestop(handle)
     182           72 :    END SUBROUTINE prepare_per_atom
     183              : 
     184              : ! **************************************************************************************************
     185              : !> \brief Build the operator for the position perturbation
     186              : !> \param dcdr_env ...
     187              : !> \param qs_env ...
     188              : !> \authors Sandra Luber
     189              : !>          Edward Ditler
     190              : !>          Ravi Kumar
     191              : !>          Rangsiman Ketkaew
     192              : ! **************************************************************************************************
     193          216 :    SUBROUTINE dcdr_build_op_dR(dcdr_env, qs_env)
     194              : 
     195              :       TYPE(dcdr_env_type)                                :: dcdr_env
     196              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     197              : 
     198              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'dcdr_build_op_dR'
     199              :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     200              : 
     201              :       INTEGER                                            :: handle, ispin, nao, nmo
     202              :       TYPE(cp_fm_type)                                   :: buf
     203          216 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: opdr_sym
     204              : 
     205          216 :       CALL timeset(routineN, handle)
     206              : 
     207          216 :       nao = dcdr_env%nao
     208              : 
     209              :       ! allocate matrix for the sum of the perturbation terms of the operator (dbcsr matrix)
     210          216 :       NULLIFY (opdr_sym)
     211          216 :       CALL dbcsr_allocate_matrix_set(opdr_sym, 1)
     212          216 :       ALLOCATE (opdr_sym(1)%matrix)
     213          216 :       CALL dbcsr_copy(opdr_sym(1)%matrix, dcdr_env%matrix_s1(1)%matrix)  ! symmetric
     214          216 :       CALL dbcsr_set(opdr_sym(1)%matrix, 0.0_dp)
     215              : 
     216          468 :       DO ispin = 1, dcdr_env%nspins
     217          252 :          nmo = dcdr_env%nmo(ispin)
     218              : 
     219          252 :          CALL apply_op_constant_term(qs_env, dcdr_env)  ! dcdr_env%matrix_apply_op_constant
     220              :          ! Hartree and Exchange-Correlation contributions
     221          252 :          CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_core_charge_1(dcdr_env%beta)%matrix, zero, one)
     222          252 :          CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_d_vhxc_dR(dcdr_env%beta, ispin)%matrix, one, one)
     223          252 :          CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_vhxc_perturbed_basis(ispin, dcdr_env%beta)%matrix, one, one)
     224              : 
     225              :          ! Core Hamiltonian contributions
     226          252 :          CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_hc(dcdr_env%beta)%matrix, one, one)
     227          252 :          CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_ppnl_1(dcdr_env%beta)%matrix, one, one)
     228          252 :          CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_apply_op_constant(ispin)%matrix, one, one)
     229              : 
     230          252 :          CALL dbcsr_desymmetrize(opdr_sym(1)%matrix, dcdr_env%hamiltonian1(1)%matrix)
     231          252 :          CALL dbcsr_add(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%matrix_t1(dcdr_env%beta + 1)%matrix, one, one)
     232              : 
     233              :          CALL cp_dbcsr_sm_fm_multiply(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%mo_coeff(ispin), &
     234          252 :                                       dcdr_env%op_dR(ispin), ncol=nmo)
     235              : 
     236              :          ! The overlap derivative terms for the Sternheimer equation
     237              :          ! buf = mo * (-mo * matrix_ks * mo)
     238          252 :          CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
     239              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
     240              :                             -1.0_dp, dcdr_env%mo_coeff(ispin), dcdr_env%chc(ispin), &
     241          252 :                             0.0_dp, buf)
     242              : 
     243              :          CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, buf, dcdr_env%op_dR(ispin), &
     244          252 :                                       nmo, alpha=1.0_dp, beta=1.0_dp)
     245          252 :          CALL cp_fm_release(buf)
     246              : 
     247              :          ! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
     248          252 :          CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(ispin))
     249              : 
     250          720 :          IF (dcdr_env%z_matrix_method) THEN
     251           54 :             CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin))
     252              :          END IF
     253              : 
     254              :       END DO
     255              : 
     256          216 :       CALL dbcsr_deallocate_matrix_set(opdr_sym)
     257              : 
     258          216 :       CALL timestop(handle)
     259          216 :    END SUBROUTINE dcdr_build_op_dR
     260              : 
     261              : ! **************************************************************************************************
     262              : !> \brief Get the dC/dR by solving the Sternheimer equation, using the op_dR matrix
     263              : !> \param dcdr_env ...
     264              : !> \param p_env ...
     265              : !> \param qs_env ...
     266              : !> \authors SL, ED
     267              : ! **************************************************************************************************
     268          180 :    SUBROUTINE dcdr_response_dR(dcdr_env, p_env, qs_env)
     269              : 
     270              :       TYPE(dcdr_env_type)                                :: dcdr_env
     271              :       TYPE(qs_p_env_type)                                :: p_env
     272              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     273              : 
     274              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'dcdr_response_dR'
     275              : 
     276              :       INTEGER                                            :: handle, ispin, output_unit
     277              :       LOGICAL                                            :: should_stop
     278          180 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: h1_psi0, psi0_order, psi1
     279              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     280              :       TYPE(cp_logger_type), POINTER                      :: logger
     281              :       TYPE(linres_control_type), POINTER                 :: linres_control
     282          180 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     283              :       TYPE(section_vals_type), POINTER                   :: lr_section
     284              : 
     285          180 :       CALL timeset(routineN, handle)
     286          180 :       NULLIFY (linres_control, lr_section, logger)
     287              : 
     288              :       CALL get_qs_env(qs_env=qs_env, &
     289              :                       linres_control=linres_control, &
     290          180 :                       mos=mos)
     291              : 
     292          180 :       logger => cp_get_default_logger()
     293          180 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     294              : 
     295              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     296          180 :                                          extension=".linresLog")
     297          180 :       IF (output_unit > 0) THEN
     298              :          WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
     299           90 :             "*** Self consistent optimization of the response wavefunction ***"
     300              :       END IF
     301              : 
     302              :       ! allocate the vectors
     303          738 :       ALLOCATE (psi0_order(dcdr_env%nspins))
     304          738 :       ALLOCATE (psi1(dcdr_env%nspins))
     305          738 :       ALLOCATE (h1_psi0(dcdr_env%nspins))
     306              : 
     307          378 :       DO ispin = 1, dcdr_env%nspins
     308          198 :          CALL cp_fm_create(psi1(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
     309          198 :          CALL cp_fm_create(h1_psi0(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
     310          198 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     311          378 :          psi0_order(ispin) = mo_coeff
     312              :       END DO
     313              : 
     314              :       ! Restart
     315          180 :       IF (linres_control%linres_restart) THEN
     316           18 :          CALL dcdr_read_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
     317              :       ELSE
     318          342 :          DO ispin = 1, dcdr_env%nspins
     319          342 :             CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
     320              :          END DO
     321              :       END IF
     322              : 
     323          180 :       IF (output_unit > 0) THEN
     324              :          WRITE (output_unit, "(T10,A,I4,A)") &
     325           90 :             "Response to the perturbation operator referring to atom ", dcdr_env%lambda, &
     326          180 :             " displaced in "//ACHAR(dcdr_env%beta + 119)
     327              :       END IF
     328          378 :       DO ispin = 1, dcdr_env%nspins
     329          198 :          CALL cp_fm_set_all(dcdr_env%dCR(ispin), 0.0_dp)
     330          378 :          CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), h1_psi0(ispin))
     331              :       END DO
     332              : 
     333          180 :       linres_control%lr_triplet = .FALSE. ! we do singlet response
     334          180 :       linres_control%do_kernel = .TRUE.
     335          180 :       linres_control%converged = .FALSE.
     336              : 
     337              :       ! Position perturbation to get dCR
     338              :       ! (H0-E0) psi1 = (H1-E1) psi0
     339              :       ! psi1 = the perturbed wavefunction
     340              :       ! h1_psi0 = (H1-E1-S1*\varepsilon)
     341              :       ! psi0_order = the unperturbed wavefunction
     342              :       CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
     343          180 :                          output_unit, should_stop)
     344          378 :       DO ispin = 1, dcdr_env%nspins
     345          378 :          CALL cp_fm_to_fm(psi1(ispin), dcdr_env%dCR(ispin))
     346              :       END DO
     347              : 
     348              :       ! Write the new result to the restart file
     349          180 :       IF (linres_control%linres_restart) THEN
     350           18 :          CALL dcdr_write_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
     351              :       END IF
     352              : 
     353              :       ! clean up
     354          378 :       DO ispin = 1, dcdr_env%nspins
     355          198 :          CALL cp_fm_release(psi1(ispin))
     356          378 :          CALL cp_fm_release(h1_psi0(ispin))
     357              :       END DO
     358          180 :       DEALLOCATE (psi1, h1_psi0, psi0_order)
     359              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     360          180 :                                         "PRINT%PROGRAM_RUN_INFO")
     361              : 
     362          180 :       CALL timestop(handle)
     363              : 
     364          540 :    END SUBROUTINE dcdr_response_dR
     365              : 
     366              : ! **************************************************************************************************
     367              : !> \brief Calculate atomic polar tensor
     368              : !> \param qs_env ...
     369              : !> \param dcdr_env ...
     370              : !> \authors Sandra Luber
     371              : !>          Edward Ditler
     372              : !>          Ravi Kumar
     373              : !>          Rangsiman Ketkaew
     374              : ! **************************************************************************************************
     375          540 :    SUBROUTINE apt_dR(qs_env, dcdr_env)
     376              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     377              :       TYPE(dcdr_env_type)                                :: dcdr_env
     378              : 
     379              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'apt_dR'
     380              : 
     381              :       INTEGER                                            :: alpha, handle, ikind, ispin, nao, nmo
     382              :       LOGICAL                                            :: ghost
     383              :       REAL(dp)                                           :: apt_basis_derivative, &
     384              :                                                             apt_coeff_derivative, charge, f_spin, &
     385              :                                                             temp1, temp2
     386          180 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: apt_el, apt_nuc
     387              :       TYPE(cp_fm_type)                                   :: overlap1_MO, tmp_fm_like_mos
     388          180 :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: dBerry_psi0, psi1_dBerry
     389              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     390          180 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     391              :       TYPE(polar_env_type), POINTER                      :: polar_env
     392          180 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     393              : 
     394          180 :       apt_basis_derivative = 0._dp
     395          180 :       apt_coeff_derivative = 0._dp
     396              : 
     397          180 :       CALL timeset(routineN, handle)
     398              : 
     399          180 :       NULLIFY (qs_kind_set, particle_set)
     400              :       CALL get_qs_env(qs_env=qs_env, &
     401              :                       qs_kind_set=qs_kind_set, &
     402          180 :                       particle_set=particle_set)
     403              : 
     404          180 :       nao = dcdr_env%nao
     405          180 :       apt_el => dcdr_env%apt_el_dcdr
     406          180 :       apt_nuc => dcdr_env%apt_nuc_dcdr
     407              : 
     408          180 :       f_spin = 2._dp/dcdr_env%nspins
     409              : 
     410          396 :       DO ispin = 1, dcdr_env%nspins
     411              :          ! Compute S^(1,R)_(ij)
     412          216 :          CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
     413          216 :          CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
     414          216 :          nmo = dcdr_env%nmo(ispin)
     415          216 :          mo_coeff => dcdr_env%mo_coeff(ispin)
     416          216 :          CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
     417          216 :          CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
     418              :          CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
     419          216 :                                       tmp_fm_like_mos, ncol=nmo)
     420              :          CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     421              :                             1.0_dp, mo_coeff, tmp_fm_like_mos, &
     422          216 :                             0.0_dp, overlap1_MO)
     423              : 
     424              :          !   C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
     425              :          !    We get the negative of the coefficients out of the linres solver
     426              :          !    And apply the constant correction due to the overlap derivative.
     427              :          CALL parallel_gemm("N", "N", nao, nmo, nmo, &
     428              :                             -0.5_dp, mo_coeff, overlap1_MO, &
     429          216 :                             -1.0_dp, dcdr_env%dCR_prime(ispin))
     430          216 :          CALL cp_fm_release(overlap1_MO)
     431              : 
     432          864 :          DO alpha = 1, 3
     433          648 :             IF (.NOT. dcdr_env%z_matrix_method) THEN
     434              : 
     435              :                ! FIRST CONTRIBUTION: dCR * moments * mo
     436          486 :                CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
     437          486 :                CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
     438          486 :                CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
     439              :                CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
     440          486 :                               -dcdr_env%ref_point(alpha), 1._dp)
     441              : 
     442              :                CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(ispin), &
     443          486 :                                             tmp_fm_like_mos, ncol=nmo)
     444          486 :                CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
     445              : 
     446          486 :                apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
     447              :                apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
     448          486 :                   = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
     449              :             ELSE
     450          162 :                CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
     451              :                CALL get_polar_env(polar_env=polar_env, psi1_dBerry=psi1_dBerry, &
     452          162 :                                   dBerry_psi0=dBerry_psi0)
     453              : 
     454              :                ! Note that here dcdr_env%dCR_prime contains only occ-occ block contribution,
     455              :                ! dcdr_env%dCR(ispin) is zero because we didn't run response calculation for dcdR.
     456              : 
     457              :                CALL cp_fm_trace(dBerry_psi0(alpha, ispin), &
     458              :                                 dcdr_env%dCR_prime(ispin), &
     459          162 :                                 temp1)
     460              : 
     461              :                CALL cp_fm_trace(dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
     462              :                                 psi1_dBerry(alpha, ispin), &
     463          162 :                                 temp2)
     464              : 
     465          162 :                apt_coeff_derivative = temp1 - temp2
     466              : 
     467              :                ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     468              :                !  - apt_coeff_derivative , here the trace is negative to compensate the
     469              :                !                          -ve sign in APTs= - 2 Z. M_alpha
     470              : 
     471          162 :                apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
     472              :                apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
     473          162 :                   = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
     474              :             END IF
     475              : 
     476              :             ! SECOND CONTRIBUTION: We assemble all combinations of r_i, d(chi)/d(idir)
     477              :             ! difdip contains derivatives with respect to atom dcdr_env%lambda
     478              :             ! difdip(alpha, beta): < a | r_alpha | db/dR_beta >
     479              :             ! Multiply by the MO coefficients
     480          648 :             CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
     481              :             CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, mo_coeff, &
     482          648 :                                          tmp_fm_like_mos, ncol=nmo)
     483          648 :             CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_basis_derivative)
     484              : 
     485              :             ! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
     486          648 :             apt_basis_derivative = -f_spin*apt_basis_derivative
     487              :             apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
     488          864 :                apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
     489              : 
     490              :          END DO ! alpha
     491              : 
     492          612 :          CALL cp_fm_release(tmp_fm_like_mos)
     493              :       END DO !ispin
     494              : 
     495              :       ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
     496          180 :       CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
     497          180 :       CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
     498          180 :       IF (.NOT. ghost) THEN
     499              :          apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
     500          180 :             apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
     501              :       END IF
     502              : 
     503              :       ! And deallocate all the things!
     504          180 :       CALL cp_fm_release(tmp_fm_like_mos)
     505          180 :       CALL cp_fm_release(overlap1_MO)
     506              : 
     507          180 :       CALL timestop(handle)
     508          180 :    END SUBROUTINE apt_dR
     509              : 
     510              : ! **************************************************************************************************
     511              : !> \brief Calculate atomic polar tensor using the localized dipole operator
     512              : !> \param qs_env ...
     513              : !> \param dcdr_env ...
     514              : !> \authors Edward Ditler
     515              : !>          Ravi Kumar
     516              : !>          Rangsiman Ketkaew
     517              : ! **************************************************************************************************
     518           36 :    SUBROUTINE apt_dR_localization(qs_env, dcdr_env)
     519              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     520              :       TYPE(dcdr_env_type)                                :: dcdr_env
     521              : 
     522              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR_localization'
     523              : 
     524              :       INTEGER                                            :: alpha, handle, i, icenter, ikind, ispin, &
     525              :                                                             map_atom, map_molecule, &
     526              :                                                             max_nbr_center, nao, natom, nmo, &
     527              :                                                             nsubset
     528              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: mapping_atom_molecule
     529           36 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: mapping_wannier_atom
     530              :       LOGICAL                                            :: ghost
     531              :       REAL(dp)                                           :: apt_basis_derivative, &
     532              :                                                             apt_coeff_derivative, charge, f_spin, &
     533              :                                                             smallest_r, this_factor, tmp_aptcontr, &
     534              :                                                             tmp_r
     535           36 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: diagonal_elements, diagonal_elements2
     536              :       REAL(dp), DIMENSION(3)                             :: distance, r_shifted
     537           36 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: apt_el, apt_nuc
     538           36 :       REAL(dp), DIMENSION(:, :, :, :), POINTER           :: apt_center, apt_subset
     539              :       TYPE(cell_type), POINTER                           :: cell
     540           36 :       TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER        :: centers_set
     541           36 :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: dBerry_psi0, psi1_dBerry
     542              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff, overlap1_MO, tmp_fm, &
     543              :                                                             tmp_fm_like_mos, tmp_fm_momo, &
     544              :                                                             tmp_fm_momo2
     545           36 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     546           36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     547              :       TYPE(polar_env_type), POINTER                      :: polar_env
     548           36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     549              : 
     550           36 :       CALL timeset(routineN, handle)
     551              : 
     552           36 :       NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
     553              : 
     554              :       CALL get_qs_env(qs_env=qs_env, &
     555              :                       qs_kind_set=qs_kind_set, &
     556              :                       particle_set=particle_set, &
     557              :                       molecule_set=molecule_set, &
     558           36 :                       cell=cell)
     559              : 
     560           36 :       nsubset = SIZE(molecule_set)
     561           36 :       natom = SIZE(particle_set)
     562           36 :       apt_el => dcdr_env%apt_el_dcdr
     563           36 :       apt_nuc => dcdr_env%apt_nuc_dcdr
     564           36 :       apt_subset => dcdr_env%apt_el_dcdr_per_subset
     565           36 :       apt_center => dcdr_env%apt_el_dcdr_per_center
     566              : 
     567              :       ! Map wannier functions to atoms
     568           36 :       IF (dcdr_env%nspins == 1) THEN
     569           36 :          max_nbr_center = dcdr_env%nbr_center(1)
     570              :       ELSE
     571            0 :          max_nbr_center = MAX(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
     572              :       END IF
     573          144 :       ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
     574          108 :       ALLOCATE (mapping_atom_molecule(natom))
     575           36 :       centers_set => dcdr_env%centers_set
     576              : 
     577           72 :       DO ispin = 1, dcdr_env%nspins
     578          180 :          DO icenter = 1, dcdr_env%nbr_center(ispin)
     579              :             ! For every center we check which atom is closest
     580              :             CALL shift_wannier_into_cell(r=centers_set(ispin)%array(1:3, icenter), &
     581              :                                          cell=cell, &
     582          144 :                                          r_shifted=r_shifted)
     583              : 
     584          144 :             smallest_r = HUGE(0._dp)
     585          612 :             DO i = 1, natom
     586          432 :                distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
     587         1728 :                tmp_r = SUM(distance**2)
     588          576 :                IF (tmp_r < smallest_r) THEN
     589          216 :                   mapping_wannier_atom(icenter, ispin) = i
     590          216 :                   smallest_r = tmp_r
     591              :                END IF
     592              :             END DO
     593              :          END DO
     594              : 
     595              :          ! Map atoms to molecules
     596           36 :          CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
     597           72 :          IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
     598           20 :             DO icenter = 1, dcdr_env%nbr_center(ispin)
     599           16 :                map_atom = mapping_wannier_atom(icenter, ispin)
     600           20 :                map_molecule = mapping_atom_molecule(map_atom)
     601              :             END DO
     602              :          END IF
     603              :       END DO !ispin
     604              : 
     605           36 :       nao = dcdr_env%nao
     606           36 :       f_spin = 2._dp/dcdr_env%nspins
     607              : 
     608           72 :       DO ispin = 1, dcdr_env%nspins
     609              :          ! Compute S^(1,R)_(ij)
     610              : 
     611           36 :          ALLOCATE (tmp_fm_like_mos)
     612           36 :          ALLOCATE (overlap1_MO)
     613           36 :          CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
     614           36 :          CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
     615           36 :          nmo = dcdr_env%nmo(ispin)
     616           36 :          mo_coeff => dcdr_env%mo_coeff(ispin)
     617           36 :          CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
     618           36 :          CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
     619              :          CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
     620           36 :                                       tmp_fm_like_mos, ncol=nmo)
     621              :          CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     622              :                             1.0_dp, mo_coeff, tmp_fm_like_mos, &
     623           36 :                             0.0_dp, overlap1_MO)
     624              : 
     625              :          !   C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
     626              :          !    We get the negative of the coefficients out of the linres solver
     627              :          !    And apply the constant correction due to the overlap derivative.
     628              :          CALL parallel_gemm("N", "N", nao, nmo, nmo, &
     629              :                             -0.5_dp, mo_coeff, overlap1_MO, &
     630           36 :                             -1.0_dp, dcdr_env%dCR_prime(ispin))
     631           36 :          CALL cp_fm_release(overlap1_MO)
     632              : 
     633          108 :          ALLOCATE (diagonal_elements(nmo))
     634           72 :          ALLOCATE (diagonal_elements2(nmo))
     635              : 
     636              :          ! Allocate temporary matrices
     637           36 :          ALLOCATE (tmp_fm)
     638           36 :          ALLOCATE (tmp_fm_momo)
     639           36 :          ALLOCATE (tmp_fm_momo2)
     640           36 :          CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
     641           36 :          CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
     642           36 :          CALL cp_fm_create(tmp_fm_momo2, dcdr_env%momo_fm_struct(ispin)%struct)
     643              : 
     644              :          ! FIRST CONTRIBUTION: dCR * moments * mo
     645           36 :          this_factor = -2._dp*f_spin
     646          144 :          DO alpha = 1, 3
     647          108 :             IF (.NOT. dcdr_env%z_matrix_method) THEN
     648              : 
     649          540 :                DO icenter = 1, dcdr_env%nbr_center(ispin)
     650          432 :                   CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
     651              :                   CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
     652          432 :                                                  ref_point=centers_set(ispin)%array(1:3, icenter))
     653              :                   CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
     654              :                                              mo_coeff=dcdr_env%dCR_prime(ispin), work=tmp_fm, nmo=nmo, &
     655              :                                              icenter=icenter, &
     656          540 :                                              res=tmp_fm_like_mos)
     657              :                END DO
     658              : 
     659              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     660              :                                   1.0_dp, mo_coeff, tmp_fm_like_mos, &
     661          108 :                                   0.0_dp, tmp_fm_momo)
     662          108 :                CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
     663              : 
     664              :             ELSE
     665            0 :                CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
     666              :                CALL get_polar_env(polar_env=polar_env, psi1_dBerry=psi1_dBerry, &
     667            0 :                                   dBerry_psi0=dBerry_psi0)
     668              : 
     669              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     670              :                                   1.0_dp, dcdr_env%dCR_prime(ispin), dBerry_psi0(alpha, ispin), &
     671            0 :                                   0.0_dp, tmp_fm_momo)
     672            0 :                CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
     673              : 
     674              :                CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     675              :                                   1.0_dp, dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
     676            0 :                                   psi1_dBerry(alpha, ispin), 0.0_dp, tmp_fm_momo2)
     677            0 :                CALL cp_fm_get_diag(tmp_fm_momo2, diagonal_elements2)
     678              : 
     679            0 :                diagonal_elements(:) = diagonal_elements(:) - diagonal_elements2(:)
     680              :             END IF
     681              : 
     682          540 :             DO icenter = 1, dcdr_env%nbr_center(ispin)
     683          432 :                map_atom = mapping_wannier_atom(icenter, ispin)
     684          432 :                map_molecule = mapping_atom_molecule(map_atom)
     685          432 :                tmp_aptcontr = this_factor*diagonal_elements(icenter)
     686              : 
     687              :                apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
     688          432 :                   = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
     689              : 
     690              :                apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
     691          540 :                   = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
     692              :             END DO
     693              : 
     694          540 :             apt_coeff_derivative = this_factor*SUM(diagonal_elements)
     695              :             apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
     696          144 :                = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
     697              :          END DO
     698              : 
     699              :          ! SECOND CONTRIBUTION: We assemble all combinations of r_i, dphi/d(idir)
     700              :          ! build part with AOs differentiated with respect to nuclear coordinates
     701              :          ! difdip contains derivatives with respect to atom dcdr_env%lambda
     702              :          ! difdip(alpha, beta): < a | r_alpha | d b/dR_beta >
     703           36 :          this_factor = -f_spin
     704          144 :          DO alpha = 1, 3
     705          540 :             DO icenter = 1, dcdr_env%nbr_center(ispin)
     706              :                ! Build the AO matrix with the right wannier center as reference point
     707          432 :                CALL dbcsr_set(dcdr_env%matrix_difdip(1, dcdr_env%beta)%matrix, 0._dp)
     708          432 :                CALL dbcsr_set(dcdr_env%matrix_difdip(2, dcdr_env%beta)%matrix, 0._dp)
     709          432 :                CALL dbcsr_set(dcdr_env%matrix_difdip(3, dcdr_env%beta)%matrix, 0._dp)
     710              :                CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, &
     711          432 :                                     1, centers_set(ispin)%array(1:3, icenter))
     712              :                CALL multiply_localization(ao_matrix=dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, &
     713              :                                           mo_coeff=mo_coeff, work=tmp_fm, nmo=nmo, &
     714              :                                           icenter=icenter, &
     715          540 :                                           res=tmp_fm_like_mos)
     716              :             END DO ! icenter
     717              : 
     718              :             CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     719              :                                1.0_dp, mo_coeff, tmp_fm_like_mos, &
     720          108 :                                0.0_dp, tmp_fm_momo)
     721          108 :             CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
     722              : 
     723          540 :             DO icenter = 1, dcdr_env%nbr_center(ispin)
     724          432 :                map_atom = mapping_wannier_atom(icenter, ispin)
     725          432 :                map_molecule = mapping_atom_molecule(map_atom)
     726          432 :                tmp_aptcontr = this_factor*diagonal_elements(icenter)
     727              : 
     728              :                apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
     729          432 :                   = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
     730              : 
     731              :                apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
     732          540 :                   = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
     733              :             END DO
     734              : 
     735              :             ! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
     736          540 :             apt_basis_derivative = this_factor*SUM(diagonal_elements)
     737              : 
     738              :             apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
     739          144 :                = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
     740              : 
     741              :          END DO  ! alpha
     742           36 :          DEALLOCATE (diagonal_elements)
     743           36 :          DEALLOCATE (diagonal_elements2)
     744              : 
     745           36 :          CALL cp_fm_release(tmp_fm)
     746           36 :          CALL cp_fm_release(tmp_fm_like_mos)
     747           36 :          CALL cp_fm_release(tmp_fm_momo)
     748           36 :          CALL cp_fm_release(tmp_fm_momo2)
     749           36 :          DEALLOCATE (overlap1_MO)
     750           36 :          DEALLOCATE (tmp_fm)
     751           36 :          DEALLOCATE (tmp_fm_like_mos)
     752           36 :          DEALLOCATE (tmp_fm_momo)
     753          144 :          DEALLOCATE (tmp_fm_momo2)
     754              :       END DO !ispin
     755              : 
     756              :       ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
     757           36 :       CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
     758           36 :       CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
     759           36 :       IF (.NOT. ghost) THEN
     760              :          apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
     761           36 :             apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
     762              : 
     763           36 :          map_molecule = mapping_atom_molecule(dcdr_env%lambda)
     764              :          apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) &
     765           36 :             = apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) + charge
     766              :       END IF
     767              : 
     768              :       ! And deallocate all the things!
     769              : 
     770           36 :       CALL timestop(handle)
     771          108 :    END SUBROUTINE apt_dR_localization
     772              : 
     773              : END MODULE qs_dcdr
     774              : 
        

Generated by: LCOV version 2.0-1