LCOV - code coverage report
Current view: top level - src - qs_vcd.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.2 % 383 380
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              : MODULE qs_vcd
       8              :    USE atomic_kind_types,               ONLY: get_atomic_kind
       9              :    USE cell_types,                      ONLY: cell_type
      10              :    USE commutator_rpnl,                 ONLY: build_com_mom_nl
      11              :    USE cp_control_types,                ONLY: dft_control_type
      12              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      13              :                                               dbcsr_copy,&
      14              :                                               dbcsr_desymmetrize,&
      15              :                                               dbcsr_set
      16              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      17              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale,&
      18              :                                               cp_fm_scale_and_add,&
      19              :                                               cp_fm_trace
      20              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      21              :                                               cp_fm_release,&
      22              :                                               cp_fm_set_all,&
      23              :                                               cp_fm_to_fm,&
      24              :                                               cp_fm_type
      25              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      26              :                                               cp_logger_type
      27              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      28              :                                               cp_print_key_unit_nr
      29              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      30              :                                               section_vals_type
      31              :    USE kinds,                           ONLY: dp
      32              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      33              :    USE particle_types,                  ONLY: particle_type
      34              :    USE qs_dcdr_ao,                      ONLY: hr_mult_by_delta_1d
      35              :    USE qs_environment_types,            ONLY: get_qs_env,&
      36              :                                               qs_environment_type
      37              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      38              :                                               qs_kind_type
      39              :    USE qs_linres_methods,               ONLY: linres_solver
      40              :    USE qs_linres_types,                 ONLY: linres_control_type,&
      41              :                                               vcd_env_type
      42              :    USE qs_mo_types,                     ONLY: mo_set_type
      43              :    USE qs_moments,                      ONLY: dipole_velocity_deriv
      44              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      45              :    USE qs_p_env_types,                  ONLY: qs_p_env_type
      46              :    USE qs_vcd_ao,                       ONLY: build_dSdV_matrix,&
      47              :                                               build_dcom_rpnl,&
      48              :                                               build_drpnl_matrix,&
      49              :                                               build_matrix_hr_rh,&
      50              :                                               hr_mult_by_delta_3d
      51              :    USE qs_vcd_utils,                    ONLY: vcd_read_restart,&
      52              :                                               vcd_write_restart
      53              : #include "./base/base_uses.f90"
      54              : 
      55              :    IMPLICIT NONE
      56              : 
      57              :    PRIVATE
      58              :    PUBLIC :: prepare_per_atom_vcd
      59              :    PUBLIC :: vcd_build_op_dV
      60              :    PUBLIC :: vcd_response_dV
      61              :    PUBLIC :: apt_dV
      62              :    PUBLIC :: aat_dV
      63              : 
      64              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd'
      65              : 
      66              :    REAL(dp), DIMENSION(3, 3, 3), PARAMETER  :: Levi_Civita = RESHAPE((/ &
      67              :                                                           0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
      68              :                                                           0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
      69              :                                                         0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
      70              :                                                                      (/3, 3, 3/))
      71              :    INTEGER, DIMENSION(3, 3), PARAMETER :: multipole_2d_to_1d = RESHAPE([4, 5, 6, 5, 7, 8, 6, 8, 9], [3, 3])
      72              : CONTAINS
      73              : 
      74              : ! **************************************************************************************************
      75              : !> \brief Compute I_{alpha beta}^lambda = d/dV^lambda_beta <m_alpha> = d/dV^lambda_beta < r x \dot{r} >
      76              : !>        The directions alpha, beta are stored in vcd_env%dcdr_env
      77              : !> \param vcd_env ...
      78              : !> \param qs_env ...
      79              : !> \author Edward Ditler
      80              : ! **************************************************************************************************
      81           18 :    SUBROUTINE aat_dV(vcd_env, qs_env)
      82              :       TYPE(vcd_env_type)                                 :: vcd_env
      83              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      84              : 
      85              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'aat_dV'
      86              :       INTEGER, PARAMETER                                 :: ispin = 1
      87              : 
      88              :       INTEGER                                            :: alpha, delta, gamma, handle, ikind, &
      89              :                                                             my_index, nao, nmo, nspins
      90              :       LOGICAL                                            :: ghost
      91              :       REAL(dp)                                           :: aat_prefactor, aat_tmp, charge, lc_tmp, &
      92              :                                                             tmp_trace
      93              :       REAL(dp), DIMENSION(3, 3)                          :: aat_tmp_33
      94              :       TYPE(cp_fm_type)                                   :: tmp_aomo
      95              :       TYPE(dft_control_type), POINTER                    :: dft_control
      96              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      97           18 :          POINTER                                         :: sab_all, sab_orb, sap_ppnl
      98           18 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      99           18 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     100              : 
     101           18 :       CALL timeset(routineN, handle)
     102              : 
     103              :       CALL get_qs_env(qs_env=qs_env, &
     104              :                       dft_control=dft_control, &
     105              :                       sap_ppnl=sap_ppnl, &
     106              :                       sab_orb=sab_orb, &
     107              :                       sab_all=sab_all, &
     108              :                       particle_set=particle_set, &
     109           18 :                       qs_kind_set=qs_kind_set)
     110              : 
     111           18 :       CALL cp_fm_create(tmp_aomo, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
     112              : 
     113           18 :       nspins = dft_control%nspins
     114           18 :       nmo = vcd_env%dcdr_env%nmo(ispin)
     115           18 :       nao = vcd_env%dcdr_env%nao
     116              :       ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), aat_atom => vcd_env%aat_atom_nvpt)
     117              : 
     118              :          ! I_{alpha beta}^lambda = 1/2c \sum_j^occ ...
     119           18 :          aat_prefactor = 1.0_dp!/(c_light_au * 2._dp)
     120           18 :          IF (nspins .EQ. 1) aat_prefactor = aat_prefactor*2.0_dp
     121              : 
     122              :          ! The non-PP part of the AAT consists of four contributions:
     123              :          !  (A1):  + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda)
     124              :          !  (A2):  - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta ∂_delta | nu > * (nu == lambda)
     125              :          !  (B):   - P^0 * ε_{alpha gamma delta} * < mu | r_gamma | nu > * (delta == beta) * (nu == lambda)
     126              :          !  (C):   + iP^1 * ε_{alpha gamma delta} * < mu | r_gamma ∂_delta | nu >
     127              : 
     128              :          ! (A1) + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda)
     129              :          ! (A2) - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta ∂_delta | nu > * (nu == lambda)
     130              :          ! Conjecture : It doesn't matter that the beta and gamma are swapped around!
     131              :          !               We define o = | ∂_delta nu >
     132              :          !                 and then < a | r_beta r_gamma | o > = < a | r_gamma r_beta | o>
     133              :          ! (A) + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda - nu == lambda)
     134              :          ! We have built the matrices - < mu | r_beta r_gamma ∂_delta | nu > in vcd_env%moments_der
     135              :          ! moments_der(1:9; 1:3) = moments_der(x, y, z, xx, xy, xz, yy, yz, zz;
     136              :          !                                     x, y, z)
     137              : 
     138           18 :          aat_tmp_33 = 0._dp
     139           72 :          DO gamma = 1, 3
     140           54 :             my_index = multipole_2d_to_1d(vcd_env%dcdr_env%beta, gamma)
     141          234 :             DO delta = 1, 3
     142              :                ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
     143              :                ! matrix_nosym_temp = - < mu | r_beta r_gamma ∂_delta | nu > * (mu - nu)
     144              :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     145          162 :                                vcd_env%moments_der_right(my_index, delta)%matrix)
     146              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     147              :                               vcd_env%moments_der_left(my_index, delta)%matrix, &
     148          162 :                               1._dp, -1._dp)
     149              : 
     150          162 :                CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     151          216 :                CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp_33(gamma, delta))
     152              :             END DO
     153              :          END DO
     154              : 
     155           72 :          DO alpha = 1, 3
     156           54 :             aat_tmp = 0._dp
     157              : 
     158              :             ! There are two remaining combinations for gamma and delta.
     159          216 :             DO gamma = 1, 3
     160          702 :                DO delta = 1, 3
     161          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     162          486 :                   IF (lc_tmp == 0._dp) CYCLE
     163              : 
     164              :                   ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
     165              :                   ! matrix_nosym_temp = - < mu | r_beta r_gamma ∂_delta | nu > * (mu - nu)
     166              :                   ! Because of the negative in moments_der, we need another negative sign here.
     167          648 :                   aat_tmp = aat_tmp + lc_tmp*aat_prefactor*aat_tmp_33(gamma, delta)
     168              :                END DO
     169              :             END DO
     170              : 
     171              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     172           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     173              :          END DO
     174              : 
     175              :          !  (B):   - P^0 * ε_{alpha gamma delta} * < mu | r_gamma | nu > * (delta == beta) * (nu == lambda)
     176              :          !      =  - P^0 * ε_{alpha gamma beta} * < mu | r_gamma | nu > * (nu == lambda)
     177              : 
     178           72 :          DO alpha = 1, 3
     179           54 :             aat_tmp = 0._dp
     180              : 
     181          216 :             DO gamma = 1, 3
     182          162 :                lc_tmp = Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta)
     183          162 :                IF (lc_tmp == 0._dp) CYCLE
     184              : 
     185              :                ! matrix_nosym_temp = < mu | r_gamma | nu > * (nu == lambda)
     186              :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(gamma)%matrix, &
     187           36 :                                        vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     188              :                CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     189           36 :                                         sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     190              : 
     191           36 :                CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     192           36 :                CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     193          216 :                aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace
     194              :             END DO
     195              : 
     196              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     197           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     198              :          END DO
     199              : 
     200              :          !  (C):   + iP^1 * ε_{alpha gamma delta} * < mu | r_gamma ∂_delta | nu >
     201           72 :          DO alpha = 1, 3
     202           54 :             aat_tmp = 0._dp
     203              : 
     204          216 :             DO gamma = 1, 3
     205          702 :                DO delta = 1, 3
     206          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     207          486 :                   IF (lc_tmp == 0._dp) CYCLE
     208              : 
     209          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%moments_der(gamma, delta)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     210          108 :                   CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), tmp_trace)
     211              : 
     212              :                   ! mo_coeff * dCV_prime = + iP1
     213              :                   ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
     214              :                   ! so we need the opposite sign.
     215          648 :                   aat_tmp = aat_tmp - 2._dp*aat_prefactor*tmp_trace*lc_tmp
     216              :                END DO
     217              :             END DO
     218              : 
     219              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     220           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     221              :          END DO
     222              : 
     223              :          ! The PP part consists of four contributions
     224              :          !  (D):  - P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma [V, r_delta] | nu > * (mu == lambda)
     225              :          !  (E):  + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
     226              :          !  (F):  - P^0 * ε_{alpha gamma delta} * < mu | r_gamma [[V, r_beta], r_delta] | nu > * (eta == lambda)
     227              :          !  (G):  - iP^1 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] | nu >
     228              : 
     229              :          !  (D):  - P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma [V, r_delta] | nu > * (mu == lambda)
     230              :          !    The negative of this is in vcd_env%matrix_r_rxvr
     231           72 :          DO alpha = 1, 3
     232              :             CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     233           54 :                             vcd_env%matrix_r_rxvr(alpha, vcd_env%dcdr_env%beta)%matrix)
     234              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     235           54 :                                      sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     236              : 
     237           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     238           54 :             CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp)
     239           54 :             aat_tmp = -aat_prefactor*aat_tmp
     240              : 
     241              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     242           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     243              :          END DO
     244              : 
     245              :          !  (E):  + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
     246              :          !    This is in vcd_env%matrix_rxvr_r
     247           72 :          DO alpha = 1, 3
     248           54 :            CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rxvr_r(alpha, vcd_env%dcdr_env%beta)%matrix)
     249              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     250           54 :                                      sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     251              : 
     252           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     253           54 :             CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp)
     254           54 :             aat_tmp = aat_prefactor*aat_tmp
     255              : 
     256              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     257           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     258              :          END DO
     259              : 
     260              :          !  (F):  - P^0 * ε_{alpha gamma delta} * < mu | r_gamma [[V, r_beta], r_delta] | nu > * (eta == lambda)
     261              :          !        + P^0 * ε_{alpha gamma delta} * < mu | [[V, r_beta], r_delta] | nu > * (eta == lambda) * R_gamma
     262              :          !    The negative is in vcd_env%matrix_r_doublecom
     263           72 :          DO alpha = 1, 3
     264              :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_r_doublecom(alpha, vcd_env%dcdr_env%beta)%matrix, &
     265           54 :                                          mo_coeff, tmp_aomo, ncol=nmo)
     266           54 :             CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp)
     267           54 :             aat_tmp = -aat_prefactor*aat_tmp
     268              : 
     269              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     270           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     271              :          END DO
     272              : 
     273              :          !  (G):   - iP^1 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] | nu >
     274           72 :          DO alpha = 1, 3
     275           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_rxrv(alpha)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     276           54 :             CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), aat_tmp)
     277              : 
     278              :             !  I can take the positive, because build_com_mom_nl computes r x [r, V]
     279           54 :             aat_tmp = 2._dp*aat_prefactor*aat_tmp
     280              : 
     281              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     282              :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     283           72 :                  + aat_tmp
     284              :          END DO
     285              : 
     286              :          ! All the reference dependent stuff
     287              :          ! (C) iP^1 * ε_{alpha gamma delta} * < mu | ∂_delta | nu > * (- R_gamma)
     288           72 :          DO alpha = 1, 3
     289           54 :             aat_tmp = 0._dp
     290              : 
     291          216 :             DO gamma = 1, 3
     292          702 :                DO delta = 1, 3
     293          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     294          486 :                   IF (lc_tmp == 0._dp) CYCLE
     295              :                   ! dipvel_ao = + < a | ∂ | b >
     296          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao(delta)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     297          108 :                   CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), tmp_trace)
     298              : 
     299              :                   ! The negative sign is due to (r - O^mag_gamma) and otherwise this is
     300              :                   !   exactly the APT dipvel(beta, delta) * (-O^mag_gamma)
     301          648 :                   aat_tmp = aat_tmp + 2._dp*aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
     302              :                END DO
     303              :             END DO
     304              : 
     305              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     306           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     307              :          END DO
     308              : 
     309              :          !  (G):  - iP^1 * ε_{alpha gamma delta} * < mu | [V, r_delta] | nu > * (- R_gamma)
     310           72 :          DO alpha = 1, 3
     311           54 :             aat_tmp = 0._dp
     312          216 :             DO gamma = 1, 3
     313          702 :                DO delta = 1, 3
     314          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     315          486 :                   IF (lc_tmp == 0._dp) CYCLE
     316              :                   ! hcom = < a | [r, V] | b > = - < a | [V, r] | b >
     317              :                   ! mo_coeff * dCV_prime = + iP1
     318          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%hcom(delta)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     319          108 :                   CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), tmp_trace)
     320              : 
     321              :                   ! This is exactly APT hcom(beta, delta)
     322          648 :                   aat_tmp = aat_tmp + 2._dp*aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
     323              :                END DO
     324              :             END DO
     325              : 
     326              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     327           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     328              :          END DO
     329              : 
     330              :          !  mag_vel, vel, mag
     331              :          ! Ai)   + ε_{alpha gamma delta} * R_beta R_gamma * < mu | ∂_delta | nu > * (mu - nu)
     332              :          ! Aii)  + ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma ∂_delta | nu > * (mu - nu)
     333              :          ! Aiii) + ε_{alpha gamma delta} * (-R_gamma) * < mu | r_beta ∂_delta | nu > * (mu - nu)
     334           72 :          DO alpha = 1, 3
     335           54 :             aat_tmp = 0._dp
     336          216 :             DO gamma = 1, 3
     337          702 :                DO delta = 1, 3
     338          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     339          486 :                   IF (lc_tmp == 0._dp) CYCLE
     340              :                   ! iii) - R_gamma * < mu | r_beta ∂_delta | nu > * (mu - nu)
     341              :                   ! mag
     342              :                   ! matrix_difdip2(beta, alpha) = - < a | r_beta | ∂_alpha b >  * (mu - nu)
     343              :                   !   so I need matrix_difdip2(beta, delta)
     344              :                   ! Only this part correspond to the APT difdip(beta, alpha)
     345              :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(vcd_env%dcdr_env%beta, delta)%matrix, mo_coeff, &
     346          108 :                                                tmp_aomo, ncol=nmo)
     347          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     348              : 
     349              :                   ! There is a negative sign here, because the routine dipole_velocity_deriv calculates
     350              :                   !   the derivatives with respect to nuclear positions and we need electronic positions
     351          108 :                   aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace*(-vcd_env%magnetic_origin_atom(gamma))
     352              : 
     353              :                   ! This part doesn't appear in the APT
     354              :                   ! ii)  - R_beta * < mu | r_gamma ∂_delta | nu > * (mu - nu)
     355              :                   ! vel
     356              :                   ! matrix_difdip2(beta, alpha) = - < a | r_beta | ∂_alpha b > * (mu - nu)
     357              :                   !   so I need matrix_difdip2(gamma, delta)
     358              :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(gamma, delta)%matrix, mo_coeff, &
     359          108 :                                                tmp_aomo, ncol=nmo)
     360          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     361              : 
     362              :                   ! There is a negative sign here, because the routine dipole_velocity_deriv calculates
     363              :                   !   the derivatives with respect to nuclear positions and we need electronic positions
     364          108 :                   aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace*(-vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta))
     365              : 
     366              :                   ! i)   + R_beta R_gamma * < mu | ∂_delta | nu > * (mu - nu)
     367              :                   ! mag_vel
     368              :                   ! dipvel_ao = + < a | ∂ | b >
     369          108 :                   CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     370          108 :                   CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
     371              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     372          108 :                                            sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     373              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
     374          108 :                                            sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     375              :                   CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     376          108 :                                  1._dp, -1._dp)
     377              : 
     378          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     379          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     380              :                   aat_tmp = aat_tmp + lc_tmp*aat_prefactor*tmp_trace* &
     381          864 :                             (vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
     382              : 
     383              :                END DO
     384              :             END DO
     385              : 
     386              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     387           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     388              :          END DO
     389              : 
     390              :          ! (B):  P^0 * ε_{alpha gamma beta} * < mu | nu > * (nu == lambda) * R_gamma
     391           72 :          DO alpha = 1, 3
     392           54 :             aat_tmp = 0._dp
     393              : 
     394          216 :             DO gamma = 1, 3
     395          162 :                lc_tmp = Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta)
     396          162 :                IF (lc_tmp == 0._dp) CYCLE
     397           36 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, 0.0_dp)
     398           36 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     399              :                CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", sab_all, &
     400           36 :                                         vcd_env%dcdr_env%lambda, direction_Or=.TRUE.)
     401           36 :                CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     402           36 :                CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     403              : 
     404              :                ! This is in total positive because we are calculating
     405              :                !  -1/2c * P * < a | b > * (delta == beta) * (nu == lambda) * (-R_gamma)
     406              :                ! The whole term corresponds to difdip_s
     407          216 :                aat_tmp = aat_tmp + lc_tmp*aat_prefactor*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
     408              :             END DO
     409              : 
     410              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     411           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     412              :          END DO
     413              : 
     414              :          !  (D):  - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta [V, r_delta] | nu > * (mu == lambda)
     415              :          !  mag, vel, mag_vel
     416              :          ! Di)   - ε_{alpha gamma delta} * (-R_gamma) * < mu | r_beta [V, r_delta] | nu > * (mu == lambda)
     417              :          ! Dii)  - ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (mu == lambda)
     418              :          ! Diii) - ε_{alpha gamma delta} * R_beta R_gamma * < mu | [V, r_delta] | nu > * (mu == lambda)
     419              : 
     420           72 :          DO alpha = 1, 3
     421           54 :             aat_tmp = 0._dp
     422           54 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
     423              : 
     424          216 :             DO gamma = 1, 3
     425          702 :                DO delta = 1, 3
     426          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     427          486 :                   IF (lc_tmp == 0._dp) CYCLE
     428              :                   ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
     429              : 
     430              :                   ! This corresponds to rcom
     431              :                   ! Di) mag
     432              :                   ! -(-R_gamma) * < mu | r_beta [V, r_delta] | nu > * (mu == lambda)
     433              :                   ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
     434              :                   !  so I need vcd_env%matrix_rrcom(delta, beta)
     435              :                   !  The multiplication with delta was not done for all directions
     436              :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     437          108 :                                   vcd_env%matrix_rrcom(delta, vcd_env%dcdr_env%beta)%matrix)
     438              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     439          108 :                                            sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     440          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     441          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     442              :                   ! The sign is positive in total, because we have the negative coordinate and the whole term was negative
     443          108 :                   aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*vcd_env%magnetic_origin_atom(gamma)
     444              : 
     445              :                   ! This doesn't appear in the APT formula
     446              :                   ! Dii) vel
     447              :                   ! -(-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (mu == lambda)
     448              :                   ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
     449              :                   !  so I need vcd_env%matrix_rrcom(delta, gamma)
     450              :                   !  The multiplication with delta was already done in SUBROUTINE apt_dV
     451          108 :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rrcom(delta, gamma)%matrix)
     452              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     453          108 :                                            sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     454          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     455          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     456          108 :                   aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)
     457              : 
     458              :                   ! Diii) mag_vel
     459              :                   !  - R_beta R_gamma * < mu | [V, r_delta] | nu >
     460              :                   ! hcom(delta) = - [V, r_delta]
     461          108 :                   CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     462              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     463          108 :                                            sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     464          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     465          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     466              :                   ! No need for a negative sign, because hcom already contains the negative sign.
     467              :                   aat_tmp = aat_tmp + &
     468              :                             aat_prefactor*tmp_trace*lc_tmp &
     469          864 :                             *(vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
     470              :                END DO
     471              :             END DO
     472              : 
     473              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     474           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     475              :          END DO
     476              : 
     477              :          !  (E):  + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
     478              :          !  mag, vel, mag_vel
     479              :          ! Ei)   + ε_{alpha gamma delta} * (-R_gamma) * < mu | [V, r_delta] r_beta | nu > * (nu == lambda)
     480              :          ! Eii)  + ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (nu == lambda)
     481              :          ! Eiii) + ε_{alpha gamma delta} * R_beta R_gamma * < mu | [V, r_delta] | nu > * (nu == lambda)
     482           72 :          DO alpha = 1, 3
     483           54 :             aat_tmp = 0._dp
     484           54 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
     485              : 
     486          216 :             DO gamma = 1, 3
     487          702 :                DO delta = 1, 3
     488          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     489          486 :                   IF (lc_tmp == 0._dp) CYCLE
     490              :                   ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
     491              :                   ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
     492              : 
     493              :                   ! This corresponds to rcom
     494              :                   ! Ei) mag
     495              :                   ! (-R_gamma) * < mu | [V, r_delta] r_beta | nu > * (nu == lambda)
     496              :                   ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
     497              :                   !  so I need vcd_env%matrix_rcomr(delta, beta)
     498              :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     499          108 :                                   vcd_env%matrix_rcomr(delta, vcd_env%dcdr_env%beta)%matrix)
     500              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     501          108 :                                            sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     502              : 
     503          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     504          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     505          108 :                   aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
     506              : 
     507              :                   ! This doesn't appear in the APT formula
     508              :                   ! E2) vel
     509              :                   ! (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (nu == lambda)
     510              :                   ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
     511              :                   !  so I need vcd_env%matrix_rrcom(delta, gamma)
     512          108 :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rrcom(delta, gamma)%matrix)
     513              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     514          108 :                                            sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     515              : 
     516          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     517          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     518          108 :                   aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta))
     519              : 
     520              :                   ! E3) mag_vel
     521              :                   ! R_beta R_gamma * < mu | [V, r_delta] | nu > * (nu == lambda)
     522          108 :                   CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     523              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     524          108 :                                            sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     525              : 
     526          108 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
     527          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     528              :                   ! There has to be a minus here, because hcom = [r, V] = - [V, r]
     529              :                   aat_tmp = aat_tmp - &
     530              :                             aat_prefactor*tmp_trace*lc_tmp* &
     531          864 :                             (vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
     532              :                END DO
     533              :             END DO
     534              : 
     535              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     536           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     537              :          END DO
     538              : 
     539              :          !  (F):  - P^0 * ε_{alpha gamma delta} * < mu | [[V, r_beta], r_delta] | nu > * (eta == lambda) * (-R_gamma)
     540              :          ! This corresponds to APT dcom
     541           72 :          DO alpha = 1, 3
     542           54 :             aat_tmp = 0._dp
     543              : 
     544          216 :             DO gamma = 1, 3
     545          702 :                DO delta = 1, 3
     546          486 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     547          486 :                   IF (lc_tmp == 0._dp) CYCLE
     548              :                   ! vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta) = - < mu | [ [V, r_beta], r_alpha ] | nu >
     549              :                   !  so I need matrix_dcom(delta, vcd_env%dcdr_env%beta)
     550              :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dcom(delta, vcd_env%dcdr_env%beta)%matrix, &
     551          108 :                                                mo_coeff, tmp_aomo, ncol=nmo)
     552          108 :                   CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
     553              :                   ! matrix_dcom has the negative sign and we include the negative sign of the coordinate
     554          648 :                   aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
     555              :                END DO
     556              :             END DO
     557              : 
     558              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     559           72 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     560              :          END DO
     561              : 
     562              :          ! Nuclear contribution
     563           18 :          CALL get_atomic_kind(particle_set(vcd_env%dcdr_env%lambda)%atomic_kind, kind_number=ikind)
     564           18 :          CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
     565           36 :          IF (.NOT. ghost) THEN
     566           72 :             DO alpha = 1, 3
     567           54 :                aat_tmp = 0._dp
     568          234 :                DO gamma = 1, 3
     569          162 :                   IF (Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta) == 0._dp) CYCLE
     570              :                   aat_tmp = aat_tmp + charge &
     571              :                             *Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta) &
     572           36 :                             *(particle_set(vcd_env%dcdr_env%lambda)%r(gamma) - vcd_env%magnetic_origin_atom(gamma))
     573              : 
     574              :                   aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     575          216 :                      = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     576              :                END DO
     577              :             END DO
     578              :          END IF
     579              :       END ASSOCIATE
     580              : 
     581           18 :       CALL cp_fm_release(tmp_aomo)
     582           18 :       CALL timestop(handle)
     583           54 :    END SUBROUTINE aat_dV
     584              : 
     585              : ! **************************************************************************************************
     586              : !> \brief Compute E_{alpha beta}^lambda = d/dV^lambda_beta <\mu_alpha> = d/dV^lambda_beta < \dot{r} >
     587              : !>        The directions alpha, beta are stored in vcd_env%dcdr_env
     588              : !> \param vcd_env ...
     589              : !> \param qs_env ...
     590              : !> \author Edward Ditler, Tomas Zimmermann
     591              : ! **************************************************************************************************
     592           18 :    SUBROUTINE apt_dV(vcd_env, qs_env)
     593              :       TYPE(vcd_env_type)                                 :: vcd_env
     594              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     595              : 
     596              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'apt_dV'
     597              :       INTEGER, PARAMETER                                 :: ispin = 1
     598              :       REAL(dp), PARAMETER                                :: f_spin = 2._dp
     599              : 
     600              :       INTEGER                                            :: alpha, handle, ikind, nao, nmo
     601              :       LOGICAL                                            :: ghost
     602              :       REAL(dp)                                           :: charge
     603              :       REAL(KIND=dp)                                      :: apt_dcom, apt_difdip, apt_dipvel, &
     604              :                                                             apt_hcom, apt_rcom
     605              :       TYPE(cp_fm_type)                                   :: buf, matrix_dSdV_mo
     606              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     607           18 :          POINTER                                         :: sab_all
     608           18 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     609           18 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     610              : 
     611           18 :       CALL timeset(routineN, handle)
     612              : 
     613              :       CALL get_qs_env(qs_env=qs_env, &
     614              :                       sab_all=sab_all, &
     615              :                       particle_set=particle_set, &
     616           18 :                       qs_kind_set=qs_kind_set)
     617              : 
     618           18 :       nmo = vcd_env%dcdr_env%nmo(ispin)
     619           18 :       nao = vcd_env%dcdr_env%nao
     620              : 
     621              :       ASSOCIATE (apt_el => vcd_env%apt_el_nvpt, &
     622              :                  apt_nuc => vcd_env%apt_nuc_nvpt, &
     623              :                  apt_total => vcd_env%apt_total_nvpt, &
     624              :                  mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), &
     625              :                  deltaR => vcd_env%dcdr_env%deltaR)
     626              : 
     627              :          ! build the full matrices
     628           18 :          CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
     629           18 :          CALL cp_fm_create(matrix_dSdV_mo, vcd_env%dcdr_env%momo_fm_struct(ispin)%struct)
     630              : 
     631              :          ! STEP 1: dCV contribution (dipvel + commutator)
     632              :          ! <mu|∂_alpha|nu> and <mu|[r_alpha, V]|nu> in AO basis
     633              :          ! We compute tr(c_1^* x ∂_munu x c_0) + tr(c_0 x ∂_munu x c_1)
     634              :          ! We compute tr(c_1^* x [,]_munu x c_0) + tr(c_0 x [,]_munu x c_1)
     635           18 :          CALL cp_fm_scale_and_add(0._dp, vcd_env%dCV_prime(ispin), -1._dp, vcd_env%dCV(ispin))
     636              : 
     637              :          ! Ref independent
     638              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
     639           18 :                                       buf, ncol=nmo)
     640              :          CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     641              :                             1.0_dp, mo_coeff, buf, &
     642           18 :                             0.0_dp, matrix_dSdV_mo)
     643              : 
     644              :          CALL parallel_gemm("N", "N", nao, nmo, nmo, &
     645              :                             -0.5_dp, mo_coeff, matrix_dSdV_mo, &
     646           18 :                             1.0_dp, vcd_env%dCV_prime(ispin))
     647              : 
     648              :          ! + i∂ - i[Vnl, r]
     649           72 :          DO alpha = 1, 3
     650           54 :             CALL cp_fm_set_all(buf, 0.0_dp)
     651              :             apt_dipvel = 0.0_dp
     652              : 
     653           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao(alpha)%matrix, mo_coeff, buf, ncol=nmo)
     654           54 :             CALL cp_fm_trace(buf, vcd_env%dCV_prime(ispin), apt_dipvel)
     655              :             ! dipvel_ao = + < a | ∂ | b >
     656              :             ! mo_coeff * dCV_prime = + iP1
     657           54 :             apt_dipvel = 2._dp*apt_dipvel
     658              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     659           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_dipvel
     660              :          END DO
     661              : 
     662           72 :          DO alpha = 1, 3
     663           54 :             CALL cp_fm_set_all(buf, 0.0_dp)
     664              :             apt_hcom = 0.0_dp
     665           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%hcom(alpha)%matrix, mo_coeff, buf, ncol=nmo)
     666           54 :             CALL cp_fm_trace(buf, vcd_env%dCV_prime(ispin), apt_hcom)
     667              : 
     668              :             ! hcom = < a | [r, V] | b > = - < a | [V, r] | b >
     669              :             ! mo_coeff * dCV_prime = + iP1
     670           54 :             apt_hcom = +2._dp*apt_hcom
     671              : 
     672              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     673           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_hcom
     674              :          END DO  !x/y/z
     675              : 
     676              :          ! STEP 2: basis function derivative contribution
     677              :       !! difdip_s
     678           18 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, 0.0_dp)
     679              :          CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s1(1)%matrix, &
     680           18 :                                  vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix)
     681              :          CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, qs_kind_set, "ORB", sab_all, &
     682           18 :                                   vcd_env%dcdr_env%lambda, direction_Or=.TRUE.)
     683              : 
     684              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
     685           18 :                                       buf, ncol=nmo, alpha=1._dp, beta=0._dp)
     686           18 :          CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
     687              : 
     688           18 :          apt_difdip = -f_spin*apt_difdip
     689              :          apt_el(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) &
     690           18 :             = apt_el(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) + apt_difdip
     691              : 
     692              :       !! difdip(j, idir) = < a | r_j | ∂_idir b >
     693              :       !! matrix_difdip2(beta, alpha) = < a | r_beta | ∂_alpha b >
     694           72 :          DO alpha = 1, 3 ! x/y/z for differentiated AO
     695              :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(vcd_env%dcdr_env%beta, alpha)%matrix, mo_coeff, &
     696           54 :                                          buf, ncol=nmo, alpha=1._dp, beta=0._dp)
     697              : 
     698           54 :             CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
     699              :             ! There is a negative sign here, because the routine dipole_velocity_deriv calculates
     700              :             !   the derivatives with respect to nuclear positions and we need electronic positions
     701           54 :             apt_difdip = -f_spin*apt_difdip
     702              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     703           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + apt_difdip
     704              : 
     705              :          END DO !alpha
     706              : 
     707              :          ! STEP 3: The terms r * [V, r]
     708              :          ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
     709              :          ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
     710           72 :          DO alpha = 1, 3 ! x/y/z for differentiated AO
     711           54 :             CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rcomr(alpha, vcd_env%dcdr_env%beta)%matrix)
     712           54 :             CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rrcom(alpha, vcd_env%dcdr_env%beta)%matrix)
     713              : 
     714              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     715           54 :                                      sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     716              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
     717           54 :                                      sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     718              : 
     719              :             CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     720           54 :                            1.0_dp, -1.0_dp)
     721              : 
     722           54 :             CALL cp_fm_set_all(buf, 0.0_dp)
     723           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     724           54 :             CALL cp_fm_trace(mo_coeff, buf, apt_rcom)
     725              : 
     726              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     727           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_rcom
     728              :          END DO !alpha
     729              : 
     730              :          ! STEP 4: pseudopotential derivative contribution
     731              :          ! vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta) = - < mu | [ [V, r_beta], r_alpha ] | nu >
     732           72 :          DO alpha = 1, 3 !x/y/z for differentiated AO
     733           54 :             CALL cp_fm_set_all(buf, 0.0_dp)
     734           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta)%matrix, mo_coeff, buf, ncol=nmo)
     735           54 :             CALL cp_fm_trace(mo_coeff, buf, apt_dcom)
     736              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     737           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_dcom
     738              :          END DO !alpha
     739              : 
     740              :          ! The reference point dependent terms:
     741              :       !! difdip_munu
     742              :          ! The additional term here is < a | db/dr(alpha)> * (delta_a - delta_b) * ref_point(beta)
     743              :          ! in qs_env%matrix_s1(2:4) there is < da/dR | b > = - < da/dr | b > = < a | db/dr >
     744           72 :          DO alpha = 1, 3
     745           54 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, 0._dp)
     746           54 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, 0._dp)
     747           54 :             CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(alpha + 1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
     748           54 :             CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
     749              : 
     750              :             ! < a | db/dr(alpha) > * R^lambda_beta * delta^lambda_nu
     751              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
     752           54 :                                      vcd_env%dcdr_env%lambda, direction_Or=.TRUE.)
     753              :             ! < a | db/dr(alpha) > * R^lambda_beta * delta^lambda_mu
     754              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
     755           54 :                                      vcd_env%dcdr_env%lambda, direction_Or=.FALSE.)
     756              : 
     757              :             ! < a | db/dr > * R^lambda_beta * ( delta^lambda_mu - delta^lambda_nu )
     758              :             CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, &
     759           54 :                            1._dp, -1._dp)
     760              : 
     761           54 :             CALL cp_fm_set_all(buf, 0.0_dp)
     762           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, mo_coeff, buf, ncol=nmo)
     763           54 :             CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
     764              : 
     765              :             ! And the whole contribution is
     766              :             ! - < a | db/dr > * (mu - nu) * ref_point
     767           54 :             apt_difdip = -apt_difdip*vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)
     768              : 
     769              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     770           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_difdip
     771              :          END DO
     772              : 
     773              :          ! And the additional factor to rcom
     774              :          ! < mu | [V, r] | nu > * R^lambda_beta * delta^lambda_mu
     775              :          ! - < mu | [V, r] | nu > * R^lambda_beta * delta^lambda_nu
     776              :          !
     777              :          !                  vcd_env%hcom(alpha) = - < mu | [V, r_alpha] | nu >
     778              :          ! particle_set(lambda)%r(vcd_env%dcdr_env%beta) = R^lambda_beta
     779              : 
     780           72 :          DO alpha = 1, 3
     781           54 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, 0._dp)
     782           54 :             CALL dbcsr_desymmetrize(vcd_env%hcom(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
     783           54 :             CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
     784              : 
     785              :             ! < mu | [V, r] | nu > * delta^lambda_nu
     786              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
     787           54 :                                      vcd_env%dcdr_env%lambda, direction_Or=.TRUE.)
     788              :             ! < mu | [V, r] | nu > * delta^lambda_mu
     789              :             CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
     790           54 :                                      vcd_env%dcdr_env%lambda, direction_Or=.FALSE.)
     791              : 
     792              :             CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, &
     793           54 :                            vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, -1._dp, +1._dp)
     794              : 
     795           54 :             CALL cp_fm_set_all(buf, 0.0_dp)
     796           54 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, mo_coeff, buf, ncol=nmo)
     797           54 :             CALL cp_fm_trace(mo_coeff, buf, apt_rcom)
     798           54 :             apt_rcom = -vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*apt_rcom
     799              : 
     800              :             apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     801           72 :                = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_rcom
     802              :          END DO
     803              : 
     804              :          ! STEP 5: nuclear contribution
     805              :          ASSOCIATE (atomic_kind => particle_set(vcd_env%dcdr_env%lambda)%atomic_kind)
     806           18 :             CALL get_atomic_kind(atomic_kind, kind_number=ikind)
     807           18 :             CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
     808           18 :             IF (.NOT. ghost) THEN
     809              :                apt_nuc(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) = &
     810           18 :                   apt_nuc(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) + charge
     811              :             END IF
     812              :          END ASSOCIATE
     813              : 
     814              :          ! STEP 6: deallocations
     815           18 :          CALL cp_fm_release(buf)
     816           72 :          CALL cp_fm_release(matrix_dSdV_mo)
     817              : 
     818              :       END ASSOCIATE
     819              : 
     820           18 :       CALL timestop(handle)
     821           18 :    END SUBROUTINE apt_dV
     822              : 
     823              : ! **************************************************************************************************
     824              : !> \brief Initialize the matrices for the NVPT calculation
     825              : !> \param vcd_env ...
     826              : !> \param qs_env ...
     827              : !> \author Edward Ditler
     828              : ! **************************************************************************************************
     829            6 :    SUBROUTINE prepare_per_atom_vcd(vcd_env, qs_env)
     830              :       TYPE(vcd_env_type)                                 :: vcd_env
     831              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     832              : 
     833              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_per_atom_vcd'
     834              : 
     835              :       INTEGER                                            :: handle, i, ispin, j
     836              :       TYPE(cell_type), POINTER                           :: cell
     837              :       TYPE(dft_control_type), POINTER                    :: dft_control
     838              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     839            6 :          POINTER                                         :: sab_all, sab_orb, sap_ppnl
     840            6 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     841            6 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     842              : 
     843            6 :       CALL timeset(routineN, handle)
     844              : 
     845              :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
     846              :                       sab_orb=sab_orb, sab_all=sab_all, sap_ppnl=sap_ppnl, &
     847            6 :                       qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
     848              : 
     849            6 :       IF (vcd_env%distributed_origin) THEN
     850            0 :          vcd_env%magnetic_origin_atom(:) = particle_set(vcd_env%dcdr_env%lambda)%r(:) - vcd_env%magnetic_origin(:)
     851            0 :          vcd_env%spatial_origin_atom = particle_set(vcd_env%dcdr_env%lambda)%r(:) - vcd_env%spatial_origin(:)
     852              :       END IF
     853              : 
     854              :       ! Reset the matrices
     855           12 :       DO ispin = 1, dft_control%nspins
     856           24 :          DO j = 1, 3
     857           18 :             CALL dbcsr_set(vcd_env%matrix_dSdV(j)%matrix, 0._dp)
     858           18 :             CALL dbcsr_set(vcd_env%matrix_drpnl(j)%matrix, 0._dp)
     859              : 
     860           78 :             DO i = 1, 3
     861           54 :                CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0.0_dp)
     862           72 :                CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0._dp)
     863              :             END DO
     864              :          END DO
     865            6 :          CALL cp_fm_set_all(vcd_env%op_dV(ispin), 0._dp)
     866           12 :          CALL dbcsr_set(vcd_env%matrix_hxc_dsdv(ispin)%matrix, 0._dp)
     867              :       END DO
     868              : 
     869              :       ! operator dV
     870              :       ! <mu|d/dV_beta [V, r_alpha]|nu>
     871              :       CALL build_dcom_rpnl(vcd_env%matrix_dcom, qs_kind_set, sab_orb, sap_ppnl, &
     872            6 :                            dft_control%qs_control%eps_ppnl, particle_set, vcd_env%dcdr_env%lambda)
     873              : 
     874              :       ! PP derivative
     875              :       CALL build_drpnl_matrix(vcd_env%matrix_drpnl, qs_kind_set, sab_all, sap_ppnl, &
     876            6 :                               dft_control%qs_control%eps_ppnl, particle_set, pseudoatom=vcd_env%dcdr_env%lambda)
     877              :       ! lin_mom
     878           24 :       DO i = 1, 3
     879           18 :          CALL dbcsr_set(vcd_env%dipvel_ao_delta(i)%matrix, 0._dp)
     880           24 :          CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dipvel_ao(i)%matrix)
     881              :       END DO
     882              : 
     883              :       CALL hr_mult_by_delta_3d(vcd_env%dipvel_ao_delta, qs_kind_set, "ORB", &
     884            6 :                                sab_all, vcd_env%dcdr_env%delta_basis_function, direction_Or=.TRUE.)
     885              : 
     886              :       ! dS/dV
     887              :       CALL build_dSdV_matrix(qs_env, vcd_env%matrix_dSdV, &
     888              :                              deltaR=vcd_env%dcdr_env%delta_basis_function, &
     889            6 :                              rcc=vcd_env%spatial_origin_atom)
     890              : 
     891              :       CALL dipole_velocity_deriv(qs_env, vcd_env%matrix_difdip2, 1, lambda=vcd_env%dcdr_env%lambda, &
     892            6 :                                  rc=[0._dp, 0._dp, 0._dp])
     893              :       ! AAT
     894              :       ! moments_throw: x, y, z, xx, xy, xz, yy, yz, zz
     895              :       ! moments_der:  (moment, xyz derivative)
     896              :       ! build_local_moments_der_matrix uses adbdr for calculating derivatives of the *primitive*
     897              :       !  on the right. So the resulting
     898              :       !  moments_der(moment, delta) = - < a | moment \partial_\delta | b >
     899           60 :       DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
     900          222 :          DO j = 1, 3
     901          162 :             CALL dbcsr_set(vcd_env%moments_der_right(i, j)%matrix, 0.0_dp)
     902          216 :             CALL dbcsr_set(vcd_env%moments_der_left(i, j)%matrix, 0.0_dp)
     903              :          END DO
     904              :       END DO
     905              : 
     906           60 :       DO i = 1, 9
     907          222 :          DO j = 1, 3 ! derivatives
     908          162 :             CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_right(i, j)%matrix) ! A2
     909          162 :             CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_left(i, j)%matrix)  ! A1
     910              : 
     911              :             !  - < mu | r_beta r_gamma ∂_delta | nu > * (mu/nu == lambda)
     912              :             CALL hr_mult_by_delta_1d(vcd_env%moments_der_right(i, j)%matrix, qs_kind_set, "ORB", &
     913          162 :                                      sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     914              :             CALL hr_mult_by_delta_1d(vcd_env%moments_der_left(i, j)%matrix, qs_kind_set, "ORB", &
     915          216 :                                      sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
     916              :          END DO
     917              :       END DO
     918              : 
     919           24 :       DO i = 1, 3
     920           78 :          DO j = 1, 3
     921           72 :             CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
     922              :          END DO
     923              :       END DO
     924              : 
     925              :       CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
     926              :                             particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
     927              :                             matrix_r_doublecom=vcd_env%matrix_r_doublecom, &
     928            6 :                             pseudoatom=vcd_env%dcdr_env%lambda)
     929              : 
     930            6 :       CALL timestop(handle)
     931              : 
     932            6 :    END SUBROUTINE prepare_per_atom_vcd
     933              : 
     934              : ! **************************************************************************************************
     935              : !> \brief What we are building here is the operator for the NVPT response:
     936              : !>     H0 * C1 - S0 * E0 * C1  = - op_dV
     937              : !>     linres_solver           = - [ H1 * C0 - S1 * C0 * E0 ]
     938              : !>   with
     939              : !>     H1 * C0 =   dH/dV * C0
     940              : !>               + i[∂]δ * C0
     941              : !>               - i S0 * C^(1,R)
     942              : !>               + i S0 * C0 * (C0 * S^(1,R) * C0)
     943              : !>               - S1 * C0 * E0
     944              : !>
     945              : !>     H1 * C0 = + i (Hr - rH) * C0                    [STEP 1]
     946              : !>               + i[∂]δ * C0                          [STEP 2]
     947              : !>               - i[V, r]δ * C0                       [STEP 3]
     948              : !>               - i S0 * C^(1,R)                      [STEP 4]
     949              : !>               - S1 * C0 * E0                        [STEP 5]
     950              : !> \param vcd_env ...
     951              : !> \param qs_env ...
     952              : !> \author Edward Ditler, Tomas Zimmermann
     953              : ! **************************************************************************************************
     954           18 :    SUBROUTINE vcd_build_op_dV(vcd_env, qs_env)
     955              :       TYPE(vcd_env_type)                                 :: vcd_env
     956              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     957              : 
     958              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vcd_build_op_dV'
     959              :       INTEGER, PARAMETER                                 :: ispin = 1
     960              : 
     961              :       INTEGER                                            :: handle, nao, nmo
     962              :       TYPE(cp_fm_type)                                   :: buf
     963              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     964           18 :          POINTER                                         :: sab_all
     965           18 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     966              : 
     967           18 :       CALL timeset(routineN, handle)
     968              : 
     969              :       CALL get_qs_env(qs_env=qs_env, &
     970              :                       sab_all=sab_all, &
     971           18 :                       qs_kind_set=qs_kind_set)
     972              : 
     973           18 :       nmo = vcd_env%dcdr_env%nmo(1)
     974           18 :       nao = vcd_env%dcdr_env%nao
     975              : 
     976           18 :       CALL build_matrix_hr_rh(vcd_env, qs_env, vcd_env%spatial_origin_atom)
     977              : 
     978              :       ! STEP 1: hr-rh
     979           18 :       CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_hr(ispin, vcd_env%dcdr_env%beta)%matrix)
     980           18 :       CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rh(ispin, vcd_env%dcdr_env%beta)%matrix)
     981              : 
     982              :       CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     983           18 :                                sab_all, vcd_env%dcdr_env%lambda, direction_or=.TRUE.)
     984              :       CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
     985           18 :                                sab_all, vcd_env%dcdr_env%lambda, direction_or=.FALSE.)
     986              :       CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     987              :                      vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     988           18 :                      1.0_dp, -1.0_dp)
     989              : 
     990              :       ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
     991              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, &
     992           18 :                                       vcd_env%op_dV(ispin), ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
     993              : 
     994              :          ! STEP 2: electronic momentum operator contribution
     995              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao_delta(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
     996              :                                       vcd_env%op_dV(ispin), &
     997           18 :                                       ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
     998              : 
     999              :          ! STEP 3: +dV_ppnl/dV, but build_drpnl_matrix gives the negative of dV_ppnl
    1000              :          ! The arguments (-1, 1) are swapped wrt to the hr-rh term, implying that
    1001              :          ! direction_Or and direction_hr do what they should.
    1002              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_drpnl(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
    1003              :                                       vcd_env%op_dV(ispin), &
    1004           18 :                                       ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
    1005              : 
    1006              :          ! STEP 4: - S0 * C^(1,R)
    1007           18 :          CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
    1008              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), &
    1009           18 :                                       vcd_env%op_dV(1), ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
    1010              : 
    1011              :          ! STEP 5: -S(1,V) * C0 * E0
    1012              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
    1013           18 :                                       buf, nmo, alpha=1.0_dp, beta=0.0_dp)
    1014              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
    1015              :                             -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &
    1016           18 :                             1.0_dp, vcd_env%op_dV(ispin))
    1017              : 
    1018           36 :          CALL cp_fm_release(buf)
    1019              :       END ASSOCIATE
    1020              : 
    1021              :       ! We have built op_dV but plug -op_dV into the linres_solver
    1022           18 :       CALL cp_fm_scale(-1.0_dp, vcd_env%op_dV(1))
    1023              : 
    1024              :       ! Revert the matrices
    1025           18 :       CALL build_matrix_hr_rh(vcd_env, qs_env, [0._dp, 0._dp, 0._dp])
    1026              : 
    1027           18 :       CALL timestop(handle)
    1028           36 :    END SUBROUTINE vcd_build_op_dV
    1029              : 
    1030              : ! *****************************************************************************
    1031              : !> \brief Get the dC/dV using the vcd_env%op_dV
    1032              : !> \param vcd_env ...
    1033              : !> \param p_env ...
    1034              : !> \param qs_env ...
    1035              : !> \author Edward Ditler, Tomas Zimmermann
    1036              : ! **************************************************************************************************
    1037           18 :    SUBROUTINE vcd_response_dV(vcd_env, p_env, qs_env)
    1038              : 
    1039              :       TYPE(vcd_env_type)                                 :: vcd_env
    1040              :       TYPE(qs_p_env_type)                                :: p_env
    1041              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1042              : 
    1043              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'vcd_response_dV'
    1044              :       INTEGER, PARAMETER                                 :: ispin = 1
    1045              : 
    1046              :       INTEGER                                            :: handle, output_unit
    1047              :       LOGICAL                                            :: failure, should_stop
    1048           72 :       TYPE(cp_fm_type), DIMENSION(1)                     :: h1_psi0, psi1
    1049              :       TYPE(cp_logger_type), POINTER                      :: logger
    1050              :       TYPE(linres_control_type), POINTER                 :: linres_control
    1051           18 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1052              :       TYPE(section_vals_type), POINTER                   :: lr_section, vcd_section
    1053              : 
    1054           18 :       CALL timeset(routineN, handle)
    1055           18 :       failure = .FALSE.
    1056              : 
    1057           18 :       NULLIFY (linres_control, lr_section, logger)
    1058              : 
    1059              :       CALL get_qs_env(qs_env=qs_env, &
    1060              :                       linres_control=linres_control, &
    1061           18 :                       mos=mos)
    1062              : 
    1063           18 :       logger => cp_get_default_logger()
    1064           18 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
    1065              :       vcd_section => section_vals_get_subs_vals(qs_env%input, &
    1066           18 :                                                 "PROPERTIES%LINRES%vcd")
    1067              : 
    1068              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
    1069           18 :                                          extension=".linresLog")
    1070           18 :       IF (output_unit > 0) THEN
    1071              :          WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
    1072            9 :             "*** Self consistent optimization of the response wavefunction ***"
    1073              :       END IF
    1074              : 
    1075              :       ASSOCIATE (psi0_order => vcd_env%dcdr_env%mo_coeff)
    1076           18 :          CALL cp_fm_create(psi1(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
    1077           18 :          CALL cp_fm_create(h1_psi0(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
    1078              : 
    1079              :          ! Restart
    1080           18 :          IF (linres_control%linres_restart) THEN
    1081           18 :             CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
    1082              :          ELSE
    1083            0 :             CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
    1084              :          END IF
    1085              : 
    1086           18 :          IF (output_unit > 0) THEN
    1087              :             WRITE (output_unit, *) &
    1088            9 :                "Response to the perturbation operator referring to the velocity of atom ", &
    1089           18 :                vcd_env%dcdr_env%lambda, " in "//ACHAR(vcd_env%dcdr_env%beta + 119)
    1090              :          END IF
    1091              : 
    1092              :          ! First response to get dCR
    1093              :          ! (H0-E0) psi1 = (H1-E1) psi0
    1094              :          ! psi1 = the perturbed wavefunction
    1095              :          ! h1_psi0 = (H1-E1)
    1096              :          ! psi0_order = the unperturbed wavefunction
    1097              :          ! Second response to get dCV
    1098           18 :          CALL cp_fm_set_all(vcd_env%dCV(ispin), 0.0_dp)
    1099           18 :          CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
    1100           18 :          CALL cp_fm_to_fm(vcd_env%op_dV(ispin), h1_psi0(ispin))
    1101              : 
    1102           18 :          linres_control%lr_triplet = .FALSE. ! we do singlet response
    1103           18 :          linres_control%do_kernel = .FALSE. ! no coupled response since imaginary perturbation
    1104           18 :          linres_control%converged = .FALSE.
    1105              :          CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
    1106           18 :                             output_unit, should_stop)
    1107           18 :          CALL cp_fm_to_fm(psi1(ispin), vcd_env%dCV(ispin))
    1108              : 
    1109              :          ! Write the new result to the restart file
    1110           36 :          IF (linres_control%linres_restart) THEN
    1111           18 :             CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
    1112              :          END IF
    1113              : 
    1114              :       END ASSOCIATE
    1115              : 
    1116              :       ! clean up
    1117           18 :       CALL cp_fm_release(psi1(ispin))
    1118           18 :       CALL cp_fm_release(h1_psi0(ispin))
    1119              : 
    1120              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
    1121           18 :                                         "PRINT%PROGRAM_RUN_INFO")
    1122              : 
    1123           18 :       CALL timestop(handle)
    1124           36 :    END SUBROUTINE vcd_response_dV
    1125              : 
    1126              : END MODULE qs_vcd
        

Generated by: LCOV version 2.0-1