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

            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_mfp
       8              :    USE atomic_kind_types,               ONLY: get_atomic_kind
       9              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      10              :                                               gto_basis_set_type
      11              :    USE cell_types,                      ONLY: cell_type
      12              :    USE commutator_rpnl,                 ONLY: build_com_vnl_giao
      13              :    USE cp_control_types,                ONLY: dft_control_type
      14              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      15              :                                               dbcsr_copy,&
      16              :                                               dbcsr_desymmetrize,&
      17              :                                               dbcsr_get_block_p,&
      18              :                                               dbcsr_p_type,&
      19              :                                               dbcsr_scale,&
      20              :                                               dbcsr_set,&
      21              :                                               dbcsr_type
      22              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      23              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_trace
      24              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25              :                                               cp_fm_release,&
      26              :                                               cp_fm_set_all,&
      27              :                                               cp_fm_to_fm,&
      28              :                                               cp_fm_type
      29              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      30              :                                               cp_logger_type
      31              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      32              :                                               cp_print_key_unit_nr
      33              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      34              :                                               section_vals_type
      35              :    USE kinds,                           ONLY: dp
      36              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      37              :    USE particle_types,                  ONLY: particle_type
      38              :    USE qs_environment_types,            ONLY: get_qs_env,&
      39              :                                               qs_environment_type
      40              :    USE qs_integral_utils,               ONLY: basis_set_list_setup
      41              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      42              :                                               qs_kind_type
      43              :    USE qs_linres_methods,               ONLY: linres_solver
      44              :    USE qs_linres_types,                 ONLY: linres_control_type,&
      45              :                                               vcd_env_type
      46              :    USE qs_mo_types,                     ONLY: mo_set_type
      47              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      48              :                                               get_neighbor_list_set_p,&
      49              :                                               neighbor_list_iterate,&
      50              :                                               neighbor_list_iterator_create,&
      51              :                                               neighbor_list_iterator_p_type,&
      52              :                                               neighbor_list_iterator_release,&
      53              :                                               neighbor_list_set_p_type
      54              :    USE qs_p_env_types,                  ONLY: qs_p_env_type
      55              :    USE qs_vcd_ao,                       ONLY: hr_mult_by_delta_1d
      56              :    USE qs_vcd_utils,                    ONLY: vcd_read_restart,&
      57              :                                               vcd_write_restart
      58              : 
      59              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      60              : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
      61              : !$                    omp_init_lock, omp_set_lock, &
      62              : !$                    omp_unset_lock, omp_destroy_lock
      63              : #include "./base/base_uses.f90"
      64              : 
      65              :    IMPLICIT NONE
      66              : 
      67              :    PRIVATE
      68              :    PUBLIC :: mfp_build_operator_gauge_independent
      69              :    PUBLIC :: mfp_build_operator_gauge_dependent
      70              :    PUBLIC :: mfp_response
      71              :    PUBLIC :: mfp_aat
      72              :    PUBLIC :: multiply_by_position
      73              : 
      74              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mfp'
      75              : 
      76              :    REAL(dp), DIMENSION(3, 3, 3), PARAMETER  :: Levi_Civita = RESHAPE((/ &
      77              :                                                           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, &
      78              :                                                           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, &
      79              :                                                         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/), &
      80              :                                                                      (/3, 3, 3/))
      81              : CONTAINS
      82              : 
      83              : ! **************************************************************************************************
      84              : !> \brief ...
      85              : !> \param vcd_env ...
      86              : !> \param qs_env ...
      87              : !> \author Edward Ditler
      88              : ! **************************************************************************************************
      89            0 :    SUBROUTINE mfp_aat(vcd_env, qs_env)
      90              :       TYPE(vcd_env_type)                                 :: vcd_env
      91              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      92              : 
      93              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'mfp_aat'
      94              :       INTEGER, PARAMETER                                 :: ispin = 1
      95              :       REAL(dp), DIMENSION(3), PARAMETER                  :: gauge_origin = 0._dp
      96              :       REAL(dp), PARAMETER                                :: f_spin = 2._dp
      97              : 
      98              :       INTEGER                                            :: alpha, delta, gamma, handle, ikind, nao, &
      99              :                                                             nmo
     100              :       LOGICAL                                            :: ghost
     101              :       REAL(dp)                                           :: aat_linmom, aat_moment, aat_moment_der, &
     102              :                                                             aat_overlap, aat_tmp, charge, lc_tmp, &
     103              :                                                             tmp_trace
     104              :       TYPE(cp_fm_type)                                   :: buf, buf2, matrix_dSdB_mo
     105              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     106            0 :          POINTER                                         :: sab_all
     107            0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     108            0 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     109              : 
     110            0 :       CALL timeset(routineN, handle)
     111              : 
     112              :       ! Some setup
     113            0 :       nmo = vcd_env%dcdr_env%nmo(ispin)
     114            0 :       nao = vcd_env%dcdr_env%nao
     115              : 
     116              :       ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), aat_atom => vcd_env%aat_atom_mfp)
     117              : 
     118              :          CALL get_qs_env(qs_env=qs_env, &
     119              :                          sab_all=sab_all, &
     120              :                          qs_kind_set=qs_kind_set, &
     121            0 :                          particle_set=particle_set)
     122              : 
     123              :          ! Set up the matrices
     124            0 :          CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
     125            0 :          CALL cp_fm_create(buf2, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
     126            0 :          CALL cp_fm_create(matrix_dSdB_mo, vcd_env%dcdr_env%momo_fm_struct(ispin)%struct)
     127              : 
     128              :          ! First prepare the coefficients dCB_prime = -dCB - 0.5 S1_ij
     129              :          !   The first negative because we get the negative coefficients out of the linres solver
     130              :          !   The second term due to the occ-occ block of the U matrix in C1 = U * C0
     131            0 :          DO alpha = 1, 3
     132              :             ! CALL cp_fm_scale_and_add(0._dp, vcd_env%dCB_prime(alpha), -1._dp, vcd_env%dCB(alpha))
     133              :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdB(alpha)%matrix, mo_coeff, &
     134            0 :                                          buf, ncol=nmo, alpha=1._dp, beta=0._dp)
     135              :             CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     136              :                                1.0_dp, mo_coeff, buf, &
     137            0 :                                0.0_dp, matrix_dSdB_mo)
     138              : 
     139              :             CALL parallel_gemm("N", "N", nao, nmo, nmo, &
     140              :                                -0.5_dp, mo_coeff, matrix_dSdB_mo, &
     141            0 :                                0.0_dp, vcd_env%dCB_prime(alpha))
     142              : 
     143              :          END DO
     144              : 
     145              :          ! The AATs in MFP consist of four terms
     146              :          !
     147              :          ! i)     (i CB_alpha) * CR_beta * S0
     148              :          ! ii)    (i CB_alpha) * C0 * < mu | d/dR_beta^\lambda nu >
     149              :          ! iii)  -1/c * ε_{alpha gamma delta} * C0 * C0 * < mu | r_delta | d/dR_beta^\lambda nu > * R^mu_gamma
     150              :          ! iv)   -1/c * ε_{alpha gamma delta} * C0 * CR_beta * < mu | r_delta | nu > * R^mu_gamma
     151              : 
     152              :          ! iii)  -1/c * ε_{alpha gamma delta} * C0 * C0 * < mu | r_delta | d/dR_beta^\lambda nu > * R^mu_gamma
     153              :          !       +1/c * ε_{alpha gamma delta} * C0 * C0 * < mu | r_delta ∂_beta | nu > * R^mu_gamma * delta_nu^lambda
     154              :          !
     155              :          !
     156              :          ! vcd_env%moments_der_right(delta, beta)%matrix
     157              :          !    = +/- < mu | r_delta ∂_beta | nu > * (nu == lambda)
     158            0 :          DO alpha = 1, 3
     159            0 :             aat_moment_der = 0._dp
     160              : 
     161            0 :             DO gamma = 1, 3
     162            0 :                DO delta = 1, 3
     163            0 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     164            0 :                   IF (lc_tmp == 0._dp) CYCLE
     165              : 
     166              :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     167            0 :                                   vcd_env%moments_der_right(delta, vcd_env%dcdr_env%beta)%matrix)
     168              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     169              :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     170            0 :                                             gauge_origin=gauge_origin)
     171              : 
     172            0 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     173            0 :                   CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
     174              : 
     175            0 :                   aat_moment_der = aat_moment_der - lc_tmp*tmp_trace
     176              :                END DO
     177              :             END DO
     178              : 
     179              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     180            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_moment_der
     181              :          END DO
     182              : 
     183              :          ! And additionally we have the reference dependence
     184              :          !  + ε_{alpha gamma delta} * P0 * < mu | R^mu_gamma * r_delta ∂_beta | nu > * delta_nu^lambda
     185              :          !  - ε_{alpha gamma delta} * P0 * < mu | r_delta ∂_beta | nu >              * R^G_gamma * delta_nu^lambda
     186              :          !  - ε_{alpha gamma delta} * P0 * < mu | R^mu_gamma * ∂_beta | nu >         * R^G_delta * delta_nu^lambda
     187              :          !  + ε_{alpha gamma delta} * P0 * < mu | ∂_beta | nu >                      * R^G_gamma * R^G_delta * delta_nu^lambda
     188            0 :          DO alpha = 1, 3
     189            0 :             aat_tmp = 0._dp
     190            0 :             DO gamma = 1, 3
     191            0 :                DO delta = 1, 3
     192            0 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     193            0 :                   IF (lc_tmp == 0._dp) CYCLE
     194              : 
     195              :                   ! - < mu | r_delta ∂_beta | nu > * R^G_gamma * delta_nu^lambda
     196              :                   ! I have
     197              :                   !  vcd_env%moments_der_right(delta, beta)%matrix = -< mu | r_delta ∂_beta | nu > * (nu == lambda)
     198              :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%moments_der_right(delta, vcd_env%dcdr_env%beta)%matrix, &
     199            0 :                                                mo_coeff, buf, ncol=nmo)
     200            0 :                   CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
     201            0 :                   aat_tmp = aat_tmp + lc_tmp*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
     202              : 
     203              :                   ! - < mu | R^mu_gamma * ∂_beta | nu > * R^G_delta * delta_nu^lambda
     204              :                   ! I have
     205              :                   !  dipvel_ao(beta) = < a | ∂_beta | b >
     206            0 :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dipvel_ao(vcd_env%dcdr_env%beta)%matrix)
     207              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     208            0 :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE.)
     209              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     210            0 :                                            sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     211              : 
     212              :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, &
     213            0 :                                                ncol=nmo, alpha=1._dp, beta=0._dp)
     214            0 :                   CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
     215            0 :                   aat_tmp = aat_tmp - lc_tmp*tmp_trace*vcd_env%spatial_origin_atom(delta)
     216              : 
     217              :                   ! + < mu | ∂_beta | nu > * R^G_gamma * R^G_delta * delta_nu^lambda
     218              :                   ! I have
     219              :                   ! dipvel_ao(beta) = < a | ∂_beta | b >
     220            0 :                   CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dipvel_ao(vcd_env%dcdr_env%beta)%matrix)
     221              :                   CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     222            0 :                                            sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     223              :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, &
     224            0 :                                                ncol=nmo, alpha=1._dp, beta=0._dp)
     225            0 :                   CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
     226              :                   aat_tmp = aat_tmp + lc_tmp*tmp_trace &
     227            0 :                             *vcd_env%spatial_origin_atom(delta)*vcd_env%magnetic_origin_atom(gamma)
     228              : 
     229              :                END DO
     230              :             END DO
     231              : 
     232              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     233            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
     234              :          END DO
     235              : 
     236              :          ! ii)    (i CB_alpha) * C0 * < mu | d/dR_beta^\lambda nu >
     237              :          !         = - (i CB_alpha) * C0 * < mu | ∂_beta | nu > * delta_nu^lambda
     238              :          !
     239              :          ! dipvel_ao = + < a | ∂ | b >, so I can take that and multiply with the delta
     240            0 :          CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(vcd_env%dcdr_env%beta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     241              :          CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
     242            0 :                                   sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
     243              : 
     244            0 :          DO alpha = 1, 3
     245            0 :             CALL cp_fm_set_all(buf, 0.0_dp)
     246              :             aat_linmom = 0.0_dp
     247              : 
     248              :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     249            0 :                                          mo_coeff, buf, ncol=nmo, alpha=1._dp, beta=0._dp)
     250            0 :             CALL cp_fm_trace(buf, vcd_env%dCB(alpha), aat_linmom)
     251              : 
     252              :             ! the minus, because dipvel_ao is the derivative wrt to the electronic coordinate,
     253              :             !   but we need the nuclear coordinate here.
     254            0 :             aat_linmom = -f_spin*aat_linmom
     255              : 
     256              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     257            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
     258              :          END DO
     259              : 
     260            0 :          DO alpha = 1, 3
     261            0 :             CALL cp_fm_set_all(buf, 0.0_dp)
     262              :             aat_linmom = 0.0_dp
     263              : 
     264            0 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     265            0 :             CALL cp_fm_trace(buf, vcd_env%dCB_prime(alpha), aat_linmom)
     266              : 
     267              :             ! the minus, because dipvel_ao is the derivative wrt to the electronic coordinate,
     268              :             !   but we need the nuclear coordinate here.
     269            0 :             aat_linmom = -f_spin*aat_linmom
     270              : 
     271              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     272            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
     273              :          END DO
     274              : 
     275              :          ! The reference dependent dSdB term is
     276              :          !   i/2c (- < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma) * R^G_delta
     277            0 :          DO alpha = 1, 3
     278            0 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
     279              : 
     280            0 :             DO gamma = 1, 3
     281            0 :                DO delta = 1, 3
     282            0 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     283            0 :                   IF (lc_tmp == 0._dp) CYCLE
     284            0 :                   CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
     285            0 :                   CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
     286            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
     287            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
     288              : 
     289              :                   ! < mu | nu > * R^mu_gamma
     290              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     291              :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     292            0 :                                             gauge_origin=gauge_origin)
     293              : 
     294              :                   ! < mu | nu > * R^nu_gamma
     295              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
     296              :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     297            0 :                                             gauge_origin=gauge_origin)
     298              : 
     299              :                   ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
     300              :                   CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     301            0 :                                  vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
     302              : 
     303              :                   ! and accumulate the results
     304              :                   CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, &
     305              :                                  vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     306            0 :                                  1._dp, vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
     307              :                END DO
     308              :             END DO
     309              : 
     310              :             ! With the reference dependent overlap derivative in vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix,
     311              :             !  we now build the corresponding coefficients as buf2 = -0.5 * C0 * S1(ref)
     312              :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
     313            0 :                                          buf, ncol=nmo, alpha=1._dp, beta=0._dp)
     314              :             CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     315              :                                1.0_dp, mo_coeff, buf, &
     316            0 :                                0.0_dp, matrix_dSdB_mo)
     317              : 
     318              :             CALL parallel_gemm("N", "N", nao, nmo, nmo, &
     319              :                                -0.5_dp, mo_coeff, matrix_dSdB_mo, &
     320            0 :                                0.0_dp, buf2)
     321              : 
     322              :             ! vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix still is < a | ∂ | b > * (nu==lambda)
     323            0 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     324            0 :             CALL cp_fm_trace(buf, buf2, aat_linmom)
     325              : 
     326            0 :             aat_linmom = -f_spin*aat_linmom
     327              : 
     328              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     329            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
     330              : 
     331              :          END DO
     332              : 
     333              :          ! iv)   -1/c * ε_{alpha gamma delta} * C0 * CR_beta * < mu | r_delta | nu > * R^mu_gamma
     334            0 :          DO alpha = 1, 3
     335            0 :             aat_moment = 0._dp
     336              : 
     337            0 :             DO gamma = 1, 3
     338            0 :                DO delta = 1, 3
     339            0 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     340            0 :                   IF (lc_tmp == 0._dp) CYCLE
     341              : 
     342            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     343              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     344              :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     345            0 :                                             gauge_origin=gauge_origin)
     346              : 
     347            0 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     348            0 :                   CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
     349              : 
     350            0 :                   aat_moment = aat_moment - lc_tmp*tmp_trace
     351              :                END DO
     352              :             END DO
     353              : 
     354              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     355            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_moment
     356              :          END DO
     357              : 
     358              :          ! And additionally we have the reference dependence
     359              :          ! - ε_{alpha gamma delta} * C0 * CR_beta * < mu | R^mu_gamma r_delta | nu >
     360              :          ! + ε_{alpha gamma delta} * C0 * CR_beta * < mu | r_delta | nu > * R^G_gamma
     361              :          ! + ε_{alpha gamma delta} * C0 * CR_beta * < mu | R^mu_gamma | nu > * R^G_delta
     362              :          ! - ε_{alpha gamma delta} * C0 * CR_beta * < mu | nu > * R^G_gamma * R^G_delta
     363            0 :          DO alpha = 1, 3
     364            0 :             aat_moment = 0._dp
     365              : 
     366            0 :             DO gamma = 1, 3
     367            0 :                DO delta = 1, 3
     368            0 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     369            0 :                   IF (lc_tmp == 0._dp) CYCLE
     370              : 
     371              :                   ! < mu | r_delta | nu > * R^G_gamma
     372            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     373            0 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     374            0 :                   CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
     375              : 
     376            0 :                   aat_moment = aat_moment + lc_tmp*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
     377              : 
     378              :                   ! < mu | R^mu_gamma | nu > * R^G_delta
     379            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     380              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     381            0 :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE.)
     382            0 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
     383            0 :                   CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
     384              : 
     385            0 :                   aat_moment = aat_moment + lc_tmp*tmp_trace*vcd_env%spatial_origin_atom(delta)
     386              : 
     387              :                   ! < mu | nu > * R^G_gamma * R^G_delta
     388            0 :                   CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, mo_coeff, buf, ncol=nmo)
     389            0 :                   CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
     390              : 
     391              :                   aat_moment = aat_moment + lc_tmp*tmp_trace &
     392            0 :                                *vcd_env%magnetic_origin_atom(gamma)*vcd_env%spatial_origin_atom(delta)
     393              :                END DO
     394              :             END DO
     395              : 
     396              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     397            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_moment
     398              :          END DO
     399              : 
     400              :          ! i)       2 * (iCB_alpha) * CR_beta * S0
     401              :          !
     402            0 :          DO alpha = 1, 3
     403            0 :             CALL cp_fm_set_all(buf, 0.0_dp)
     404              :             aat_overlap = 0.0_dp
     405              : 
     406            0 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), buf, ncol=nmo)
     407            0 :             CALL cp_fm_trace(buf, vcd_env%dCB(alpha), aat_overlap)
     408              : 
     409            0 :             aat_overlap = f_spin*aat_overlap
     410              : 
     411              :             ! dCB_prime * dCR_prime = + iP1
     412              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     413            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
     414              :          END DO
     415              : 
     416            0 :          DO alpha = 1, 3
     417            0 :             CALL cp_fm_set_all(buf, 0.0_dp)
     418              :             aat_overlap = 0.0_dp
     419              : 
     420            0 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), buf, ncol=nmo)
     421            0 :             CALL cp_fm_trace(buf, vcd_env%dCB_prime(alpha), aat_overlap)
     422              : 
     423            0 :             aat_overlap = f_spin*aat_overlap
     424              : 
     425              :             ! dCB_prime * dCR_prime = + iP1
     426              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     427            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
     428              :          END DO
     429              : 
     430              :          ! The reference dependent dSdB term is
     431              :          !   i/2c (- < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma) * R^G_delta
     432            0 :          DO alpha = 1, 3
     433            0 :             CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
     434              : 
     435            0 :             DO gamma = 1, 3
     436            0 :                DO delta = 1, 3
     437            0 :                   lc_tmp = Levi_Civita(alpha, gamma, delta)
     438            0 :                   IF (lc_tmp == 0._dp) CYCLE
     439            0 :                   CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
     440            0 :                   CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
     441            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
     442            0 :                   CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
     443              : 
     444              :                   ! < mu | nu > * R^mu_gamma
     445              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     446              :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     447            0 :                                             gauge_origin=gauge_origin)
     448              : 
     449              :                   ! < mu | nu > * R^nu_gamma
     450              :                   CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
     451              :                                             qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     452            0 :                                             gauge_origin=gauge_origin)
     453              : 
     454              :                   ! !! A = alpha*A + beta*B or
     455              :                   ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
     456              :                   CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     457            0 :                                  vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
     458              : 
     459              :                   ! and accumulate the results
     460              :                   CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     461            0 :                                  1._dp, -vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
     462              :                END DO
     463              :             END DO
     464              : 
     465              :             ! With the reference dependent overlap derivative in vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix,
     466              :             !  we now build the corresponding coefficients as buf2 = -0.5 * C0 * S1(ref)
     467              :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
     468            0 :                                          buf, ncol=nmo, alpha=1._dp, beta=0._dp)
     469              :             CALL parallel_gemm("T", "N", nmo, nmo, nao, &
     470              :                                1.0_dp, mo_coeff, buf, &
     471            0 :                                0.0_dp, matrix_dSdB_mo)
     472              : 
     473              :             CALL parallel_gemm("N", "N", nao, nmo, nmo, &
     474              :                                -0.5_dp, mo_coeff, matrix_dSdB_mo, &
     475            0 :                                0.0_dp, buf2)
     476              : 
     477            0 :             CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), buf, ncol=nmo)
     478            0 :             CALL cp_fm_trace(buf, buf2, aat_overlap)
     479              : 
     480            0 :             aat_overlap = f_spin*aat_overlap
     481              : 
     482              :             ! dCB_prime * dCR_prime = + iP1
     483              :             aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     484            0 :                = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
     485              :          END DO
     486              : 
     487              :          ! Nuclear contribution
     488            0 :          CALL get_atomic_kind(particle_set(vcd_env%dcdr_env%lambda)%atomic_kind, kind_number=ikind)
     489            0 :          CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
     490            0 :          IF (.NOT. ghost) THEN
     491            0 :             DO alpha = 1, 3
     492            0 :                aat_tmp = 0._dp
     493            0 :                DO gamma = 1, 3
     494            0 :                   IF (Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta) == 0._dp) CYCLE
     495              :                   aat_tmp = aat_tmp + charge &
     496              :                             *Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta) &
     497            0 :                             *(particle_set(vcd_env%dcdr_env%lambda)%r(gamma) - vcd_env%magnetic_origin_atom(gamma))
     498              : 
     499              :                   aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
     500            0 :                      = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp/4.
     501              :                END DO
     502              :             END DO
     503              :          END IF
     504              : 
     505              :       END ASSOCIATE
     506              : 
     507            0 :       CALL cp_fm_release(buf)
     508            0 :       CALL cp_fm_release(buf2)
     509            0 :       CALL cp_fm_release(matrix_dSdB_mo)
     510              : 
     511            0 :       CALL timestop(handle)
     512            0 :    END SUBROUTINE mfp_aat
     513              : 
     514              : ! **************************************************************************************************
     515              : !> \brief ...
     516              : !> \param vcd_env ...
     517              : !> \param qs_env ...
     518              : !> \param alpha ...
     519              : !> \author Edward Ditler
     520              : ! **************************************************************************************************
     521            0 :    SUBROUTINE mfp_build_operator_gauge_dependent(vcd_env, qs_env, alpha)
     522              :       TYPE(vcd_env_type)                                 :: vcd_env
     523              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     524              :       INTEGER                                            :: alpha
     525              : 
     526              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'mfp_build_operator_gauge_dependent'
     527              :       INTEGER, PARAMETER                                 :: ispin = 1
     528              : 
     529              :       INTEGER                                            :: delta, gamma, handle, nao, nmo
     530              :       REAL(dp)                                           :: eps_ppnl, lc_tmp
     531              :       REAL(dp), DIMENSION(3)                             :: gauge_origin
     532              :       TYPE(cell_type), POINTER                           :: cell
     533              :       TYPE(cp_fm_type)                                   :: buf
     534            0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     535              :       TYPE(dft_control_type), POINTER                    :: dft_control
     536              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     537            0 :          POINTER                                         :: sab_all, sap_ppnl
     538            0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     539            0 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     540              : 
     541              :       ! --------------------------------------------------------------- !
     542              :       ! The operator consists of four parts:
     543              :       !  1. R^mu x r * H^KS
     544              :       !  2. H^KS * R^nu x r
     545              :       !  3. [Vnl, (R^mu - Omag) x r]
     546              :       !  4. (r - Omag) x ∂
     547              :       ! and additionally the overlap derivative comes in as
     548              :       !  5. -dSdB * C0 * CHC
     549              :       ! --------------------------------------------------------------- !
     550              : 
     551            0 :       CALL timeset(routineN, handle)
     552              : 
     553              :       ! Some setup
     554            0 :       nmo = vcd_env%dcdr_env%nmo(ispin)
     555            0 :       nao = vcd_env%dcdr_env%nao
     556              :       ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
     557              : 
     558              :          CALL get_qs_env(qs_env=qs_env, &
     559              :                          sab_all=sab_all, &
     560              :                          qs_kind_set=qs_kind_set, &
     561              :                          particle_set=particle_set, &
     562              :                          sap_ppnl=sap_ppnl, &
     563              :                          cell=cell, &
     564              :                          dft_control=dft_control, &     ! For the eps_ppnl
     565            0 :                          matrix_ks=matrix_ks)
     566              : 
     567            0 :          gauge_origin(:) = vcd_env%magnetic_origin_atom
     568              : 
     569            0 :          eps_ppnl = dft_control%qs_control%eps_ppnl
     570              : 
     571            0 :          CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
     572            0 :          CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
     573            0 :          CALL cp_fm_set_all(buf, 0._dp)
     574              : 
     575              :          ! The first two terms in the operator are
     576              :          !    i/2c sum_{gamma delta} ε_{alpha gamma delta}
     577              :          !       (R^mu_gamma - O^mag_gamma) * < (r_delta - O^spatial) * H >
     578              :          !     - (R^nu_gamma - O^mag_gamma) * < H * (r_delta - O^spatial) >)
     579              : 
     580              :          ! they expand into
     581              :          ! i)     R^mu_gamma * matrix_rh(delta) - R^nu_gamma * matrix_hr(delta)
     582              :          ! ii)  - O^mag_gamma * (matrix_rh(delta) - matrix_hr(delta))
     583              :          ! iii) - (R^mu_gamma - R^nu_gamma) * O^spatial_delta * matrix_ks
     584              : 
     585            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     586            0 :          DO gamma = 1, 3
     587            0 :             DO delta = 1, 3
     588            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     589            0 :                IF (lc_tmp == 0._dp) CYCLE
     590              : 
     591              :                ! i)
     592              :                ! R^mu_gamma * matrix_rh(delta) - R^nu_gamma * matrix_hr(delta)
     593              :                ! R^mu_gamma * matrix_rh(delta)
     594            0 :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
     595              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     596              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     597            0 :                                          gauge_origin=gauge_origin)
     598              : 
     599              :                ! - R^nu_gamma * matrix_hr(delta)
     600            0 :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_hr(ispin, delta)%matrix)
     601              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     602              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     603            0 :                                          gauge_origin=gauge_origin)
     604              : 
     605              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     606            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
     607              : 
     608              :                ! and accumulate the results
     609              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     610            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
     611              : 
     612              :                ! ii)
     613              :                ! - O^mag_gamma * (matrix_rh(delta) - matrix_hr(delta))
     614            0 :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
     615            0 :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_hr(ispin, delta)%matrix)
     616              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     617            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
     618            0 :                CALL dbcsr_scale(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, -vcd_env%magnetic_origin_atom(gamma))
     619              : 
     620              :                ! and accumulate the results
     621              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     622            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
     623              : 
     624              :                ! iii)
     625              :                ! - (R^mu_gamma - R^nu_gamma) * O^spatial_delta * matrix_ks
     626            0 :                CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     627            0 :                CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
     628              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     629              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     630            0 :                                          gauge_origin=gauge_origin)
     631              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     632              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     633            0 :                                          gauge_origin=gauge_origin)
     634              : 
     635              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     636            0 :                               1._dp, -1._dp)
     637            0 :                CALL dbcsr_scale(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, -vcd_env%spatial_origin_atom(delta))
     638              :                ! and accumulate the results
     639              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     640            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
     641              : 
     642              :             END DO
     643              :          END DO
     644              : 
     645              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
     646            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp, beta=1.0_dp)
     647              : 
     648              :          ! The third term in the operator is
     649              :          !    i/2c sum_{gamma delta} ε_{alpha gamma delta}
     650              :          !       sum_\eta  < (R^eta_gamma - O^mag) * [Vnl^eta, r_delta] >
     651              :          !
     652              :          ! build_com_vnl_giao calculates
     653              :          !   matrix_rv(gamma, delta) =
     654              :          !       sum_\eta  < mu | R^eta_gamma * Vnl^eta * r_delta | nu >
     655              :          ! or the other way around: r_delta * Vnl^eta
     656            0 :          DO gamma = 1, 3
     657            0 :             DO delta = 1, 3
     658            0 :                CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, 0._dp)
     659            0 :                CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, 0._dp)
     660              :             END DO
     661              :          END DO
     662              :          ! V * r_delta * R^eta_gamma
     663              :          CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
     664              :                                  eps_ppnl=dft_control%qs_control%eps_ppnl, &
     665              :                                  particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp_33, &
     666            0 :                                  ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.TRUE.)
     667              : 
     668              :          ! r_delta * V * R^eta_gamma
     669              :          CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
     670              :                                  eps_ppnl=dft_control%qs_control%eps_ppnl, &
     671              :                                  particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp2_33, &
     672            0 :                                  ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.FALSE.)
     673              : 
     674            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     675            0 :          DO gamma = 1, 3
     676            0 :             DO delta = 1, 3
     677            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     678            0 :                IF (lc_tmp == 0._dp) CYCLE
     679              : 
     680              :                ! + V * r_delta * R^eta_gamma
     681              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, &
     682            0 :                               1._dp, lc_tmp)
     683              : 
     684              :                ! - r_delta * V * R^eta_gamma
     685              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, &
     686            0 :                               1._dp, -lc_tmp)
     687              :             END DO
     688              :          END DO
     689              : 
     690              :          ! Same factor 1/2 as above for Hr - rH
     691              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
     692            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp/2._dp, beta=1.0_dp)
     693              : 
     694              :          ! And the O^mag dependent term:
     695              :          !    - sum_\eta  < O^mag * [Vnl^eta, r_delta] >
     696            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     697            0 :          DO gamma = 1, 3
     698            0 :             DO delta = 1, 3
     699            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     700            0 :                IF (lc_tmp == 0._dp) CYCLE
     701              : 
     702              :                ! vcd_env%hcom(delta) = - < mu | [V, r_delta] | nu >
     703            0 :                CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     704              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     705            0 :                               1._dp, vcd_env%magnetic_origin_atom(gamma)*lc_tmp/(2._dp))
     706              :             END DO
     707              :          END DO
     708              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
     709            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
     710              : 
     711              :          ! The fourth term in the operator is
     712              :          !  -i/2c sum_{gamma delta} ε_{alpha gamma delta}
     713              :          !      (r_gamma - O^mag_gamma) * ∂_delta
     714              : 
     715              :          ! It's also the P^1 term of the NVPT AAT with the reversed sign
     716              :          !  So in the first term alpha=+1
     717            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     718            0 :          DO gamma = 1, 3
     719            0 :             DO delta = 1, 3
     720            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     721            0 :                IF (lc_tmp == 0._dp) CYCLE
     722              : 
     723              :                ! moments_der(gamma, delta) = - < a | r_gamma \partial_delta | b >
     724              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%moments_der(gamma, delta)%matrix, &
     725            0 :                               1._dp, lc_tmp/(2._dp))
     726              : 
     727              :             END DO
     728              :          END DO
     729              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
     730            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
     731              : 
     732              :          ! And the O^mag dependent term
     733              :          !   + i/2c sum_{gamma delta} ε_{alpha gamma delta}
     734              :          !      O^mag_gamma * ∂_delta
     735            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     736            0 :          DO gamma = 1, 3
     737            0 :             DO delta = 1, 3
     738            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     739            0 :                IF (lc_tmp == 0._dp) CYCLE
     740              : 
     741              :                ! dipvel_ao(delta) = + < a | ∂_delta | b >
     742            0 :                CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     743              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     744            0 :                               1._dp, vcd_env%magnetic_origin_atom(gamma)*lc_tmp/(2._dp))
     745              :             END DO
     746              :          END DO
     747              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
     748            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
     749              : 
     750              :          ! And the overlap derivative
     751              :          !   i/2c * sum_{gamma delta} ε_{alpha gamma delta}
     752              :          !     < r_delta > * (R^mu_gamma - R^nu_gamma)
     753              :          ! The reference independent part
     754            0 :          CALL dbcsr_set(vcd_env%matrix_dSdB(alpha)%matrix, 0._dp)
     755              : 
     756            0 :          DO gamma = 1, 3
     757            0 :             DO delta = 1, 3
     758            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     759            0 :                IF (lc_tmp == 0._dp) CYCLE
     760            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
     761            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 0._dp)
     762            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     763            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
     764              : 
     765              :                ! moments * R^mu_gamma
     766              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     767              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     768            0 :                                          gauge_origin=gauge_origin)
     769              : 
     770              :                ! moments * R^nu_gamma
     771              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     772              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     773            0 :                                          gauge_origin=gauge_origin)
     774              : 
     775              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     776            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
     777              : 
     778              :                ! and accumulate the results
     779              :                CALL dbcsr_add(vcd_env%matrix_dSdB(alpha)%matrix, &
     780            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
     781              :             END DO
     782              :          END DO
     783              : 
     784              :          ! the overlap derivative goes in as -S(1,B) * C0 * E0
     785              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdB(alpha)%matrix, mo_coeff, &
     786            0 :                                       buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
     787              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
     788              :                             -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &   ! the negative sign is here
     789            0 :                             1.0_dp, vcd_env%op_dB(ispin))
     790              : 
     791              :          ! And the reference dependent part of the overlap
     792            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
     793            0 :          DO gamma = 1, 3
     794            0 :             DO delta = 1, 3
     795            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     796            0 :                IF (lc_tmp == 0._dp) CYCLE
     797            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
     798            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
     799            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
     800            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
     801              : 
     802              :                ! < mu | nu > * R^mu_gamma
     803              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     804              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     805            0 :                                          gauge_origin=gauge_origin)
     806              : 
     807              :                ! < mu | nu > * R^nu_gamma
     808              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
     809              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     810            0 :                                          gauge_origin=gauge_origin)
     811              : 
     812              :                ! !! A = alpha*A + beta*B or
     813              :                ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
     814              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     815            0 :                               vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
     816              : 
     817              :                ! and accumulate the results
     818              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
     819            0 :                               1._dp, vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
     820              :             END DO
     821              :          END DO
     822              : 
     823              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
     824            0 :                                       buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
     825              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
     826              :                             -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &   ! again, the negative sign is here
     827            0 :                             1.0_dp, vcd_env%op_dB(ispin))
     828              : 
     829              :          ! CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
     830            0 :          CALL cp_fm_release(buf)
     831              : 
     832              :       END ASSOCIATE
     833              : 
     834              :       ! We have built op_dB but plug -op_dB into the linres_solver
     835              :       ! CALL cp_fm_scale(-1.0_dp, vcd_env%op_dB(ispin))
     836              : 
     837            0 :       CALL timestop(handle)
     838            0 :    END SUBROUTINE mfp_build_operator_gauge_dependent
     839              : 
     840              : ! **************************************************************************************************
     841              : !> \brief ...
     842              : !> \param vcd_env ...
     843              : !> \param qs_env ...
     844              : !> \param alpha ...
     845              : !> \author Edward Ditler
     846              : ! **************************************************************************************************
     847            0 :    SUBROUTINE mfp_build_operator_gauge_independent(vcd_env, qs_env, alpha)
     848              :       TYPE(vcd_env_type)                                 :: vcd_env
     849              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     850              :       INTEGER                                            :: alpha
     851              : 
     852              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'mfp_build_operator_gauge_independent'
     853              :       INTEGER, PARAMETER                                 :: ispin = 1
     854              : 
     855              :       INTEGER                                            :: delta, gamma, handle, nao, nmo
     856              :       REAL(dp)                                           :: eps_ppnl, lc_tmp
     857              :       REAL(dp), DIMENSION(3)                             :: gauge_origin
     858              :       TYPE(cell_type), POINTER                           :: cell
     859              :       TYPE(cp_fm_type)                                   :: buf
     860            0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     861              :       TYPE(dft_control_type), POINTER                    :: dft_control
     862              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     863            0 :          POINTER                                         :: sab_all, sap_ppnl
     864            0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     865            0 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     866              : 
     867              :       ! --------------------------------------------------------------- !
     868              :       ! The operator consists of three parts:
     869              :       !  1. (R^mu - R^nu) x r * H_KS
     870              :       !  2. [Vnl, (R^mu - R^nu) x r]
     871              :       !  3. (r - R^nu) x ∂
     872              :       ! and additionally the overlap derivative comes in as
     873              :       !  4. -dSdB * C0 * CHC
     874              :       ! --------------------------------------------------------------- !
     875              : 
     876            0 :       CALL timeset(routineN, handle)
     877              : 
     878              :       ! Some setup
     879            0 :       nmo = vcd_env%dcdr_env%nmo(ispin)
     880            0 :       nao = vcd_env%dcdr_env%nao
     881              :       ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
     882              : 
     883              :          CALL get_qs_env(qs_env=qs_env, &
     884              :                          sab_all=sab_all, &
     885              :                          qs_kind_set=qs_kind_set, &
     886              :                          particle_set=particle_set, &
     887              :                          sap_ppnl=sap_ppnl, &
     888              :                          cell=cell, &
     889              :                          dft_control=dft_control, &     ! For the eps_ppnl
     890            0 :                          matrix_ks=matrix_ks)
     891              : 
     892            0 :          gauge_origin(:) = 0._dp
     893              : 
     894            0 :          eps_ppnl = dft_control%qs_control%eps_ppnl
     895              : 
     896            0 :          CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
     897            0 :          CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
     898            0 :          CALL cp_fm_set_all(buf, 0._dp)
     899              : 
     900              :          ! The first term in the operator is
     901              :          !    i/2c sum_{gamma delta} ε_{alpha gamma delta}
     902              :          !       (R^mu_gamma - R^nu_gamma) < (r_delta - Ospatial) H >
     903              : 
     904              :          ! In vcd_second_prepare we build
     905              :          !       vcd_env%matrix_hr(1, delta) = < mu | H * r_delta | nu >
     906              :          ! and   vcd_env%matrix_rh(1, delta) = < mu | r_delta * H | nu >
     907              : 
     908              :          ! So we need
     909              :          !   (R^mu_gamma - R^nu_gamma) * vcd_env%matrix_rh
     910            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     911            0 :          DO gamma = 1, 3
     912            0 :             DO delta = 1, 3
     913            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     914            0 :                IF (lc_tmp == 0._dp) CYCLE
     915              : 
     916              :                ! The reference independent part
     917              :                ! vcd_env%matrix_rh * R^mu_gamma
     918            0 :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
     919              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     920              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     921            0 :                                          gauge_origin=gauge_origin)
     922              : 
     923              :                ! - vcd_env%matrix_rh * R^nu_gamma
     924            0 :                CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
     925              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     926              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     927            0 :                                          gauge_origin=gauge_origin)
     928              : 
     929              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     930            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
     931              : 
     932              :                ! and accumulate the results
     933              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     934            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
     935              : 
     936              :                ! And the Ospatial dependent part
     937              :                ! - (R^mu_gamma - R^nu_gamma) * Ospatial * matrix_ks
     938            0 :                CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
     939            0 :                CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
     940              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
     941              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
     942            0 :                                          gauge_origin=gauge_origin)
     943              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     944              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
     945            0 :                                          gauge_origin=gauge_origin)
     946              : 
     947              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
     948            0 :                               1._dp, -1._dp)
     949            0 :                CALL dbcsr_scale(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, -vcd_env%spatial_origin_atom(delta))
     950              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     951            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
     952              : 
     953              :             END DO
     954              :          END DO
     955              : 
     956              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
     957            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp, beta=1.0_dp)
     958              : 
     959              :          ! The second term in the operator is
     960              :          !    i/2c sum_{gamma delta} ε_{alpha gamma delta}
     961              :          !       sum_\eta  < (R^eta_gamma - R^nu_gamma) * [Vnl^eta, r_delta] >
     962              :          !
     963              :          ! build_com_vnl_giao calculates
     964              :          !   matrix_rv(gamma, delta) =
     965              :          !       sum_\eta  < mu | R^eta_gamma * Vnl^eta * r_delta | nu >
     966              :          ! or the other way around: r_delta * Vnl^eta
     967            0 :          DO gamma = 1, 3
     968            0 :             DO delta = 1, 3
     969            0 :                CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, 0._dp)
     970            0 :                CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, 0._dp)
     971              :             END DO
     972              :          END DO
     973              :          ! V * r_delta * R^eta_gamma
     974              :          CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
     975              :                                  eps_ppnl=dft_control%qs_control%eps_ppnl, &
     976              :                                  particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp_33, &
     977            0 :                                  ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.TRUE.)
     978              : 
     979              :          ! r_delta * V * R^eta_gamma
     980              :          CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
     981              :                                  eps_ppnl=dft_control%qs_control%eps_ppnl, &
     982              :                                  particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp2_33, &
     983            0 :                                  ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.FALSE.)
     984              : 
     985            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
     986            0 :          DO gamma = 1, 3
     987            0 :             DO delta = 1, 3
     988            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
     989            0 :                IF (lc_tmp == 0._dp) CYCLE
     990              : 
     991              :                ! + V * r_delta * R^eta_gamma
     992              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     993              :                               vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, &
     994            0 :                               1._dp, lc_tmp)
     995              : 
     996              :                ! - r_delta * V * R^eta_gamma
     997              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
     998              :                               vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, &
     999            0 :                               1._dp, -lc_tmp)
    1000              :             END DO
    1001              :          END DO
    1002              : 
    1003              :          ! Same factor 1/2 as above for Hr - rH
    1004              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
    1005            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp/2._dp, beta=1.0_dp)
    1006              : 
    1007              :          ! The third term in the operator is
    1008              :          !  -i/2c sum_{gamma delta} ε_{alpha gamma delta}
    1009              :          !      (r_gamma - R^nu_gamma) * ∂_delta
    1010              : 
    1011              :          ! It's also the P^1 term of the NVPT AAT with the reversed sign
    1012              :          !  So in the first term alpha=+1
    1013            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
    1014            0 :          DO gamma = 1, 3
    1015            0 :             DO delta = 1, 3
    1016            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
    1017            0 :                IF (lc_tmp == 0._dp) CYCLE
    1018              : 
    1019              :                ! moments_der(gamma, delta) = - < a | r_gamma \partial_delta | b >
    1020              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%moments_der(gamma, delta)%matrix, &
    1021            0 :                               1._dp, lc_tmp/(2._dp))
    1022              : 
    1023              :             END DO
    1024              :          END DO
    1025              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
    1026            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
    1027              : 
    1028              :          ! And the R^nu dependent term
    1029              :          !   + i/2c sum_{gamma delta} ε_{alpha gamma delta}
    1030              :          !      R^nu_gamma * ∂_delta
    1031            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
    1032            0 :          DO gamma = 1, 3
    1033            0 :             DO delta = 1, 3
    1034            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
    1035            0 :                IF (lc_tmp == 0._dp) CYCLE
    1036              : 
    1037              :                ! dipvel_ao(delta) = + < a | ∂_delta | b >
    1038            0 :                CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
    1039              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
    1040              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
    1041            0 :                                          gauge_origin=gauge_origin)
    1042              : 
    1043              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
    1044            0 :                               1._dp, lc_tmp/(2._dp))
    1045              :             END DO
    1046              :          END DO
    1047              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
    1048            0 :                                       vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
    1049              : 
    1050              :          ! And the overlap derivative
    1051              :          !   i/2c * sum_{gamma delta} ε_{alpha gamma delta}
    1052              :          !     < r_delta > * (R^mu_gamma - R^nu_gamma)
    1053              :          ! The reference independent part
    1054            0 :          CALL dbcsr_set(vcd_env%matrix_dSdB(alpha)%matrix, 0._dp)
    1055              : 
    1056            0 :          DO gamma = 1, 3
    1057            0 :             DO delta = 1, 3
    1058            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
    1059            0 :                IF (lc_tmp == 0._dp) CYCLE
    1060            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
    1061            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 0._dp)
    1062            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
    1063            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
    1064              : 
    1065              :                ! moments * R^mu_gamma
    1066              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
    1067              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
    1068            0 :                                          gauge_origin=gauge_origin)
    1069              : 
    1070              :                ! moments * R^nu_gamma
    1071              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
    1072              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
    1073            0 :                                          gauge_origin=gauge_origin)
    1074              : 
    1075              :                ! !! A = alpha*A + beta*B or
    1076              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
    1077            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
    1078              : 
    1079              :                ! and accumulate the results
    1080              :                CALL dbcsr_add(vcd_env%matrix_dSdB(alpha)%matrix, &
    1081            0 :                               vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
    1082              :             END DO
    1083              :          END DO
    1084              : 
    1085              :          ! the overlap derivative goes in as -S(1,B) * C0 * E0
    1086              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdB(alpha)%matrix, mo_coeff, &
    1087            0 :                                       buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
    1088              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
    1089              :                             -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &   ! the negative sign is here
    1090            0 :                             1.0_dp, vcd_env%op_dB(ispin))
    1091              : 
    1092              :          ! And the reference dependent part of the overlap
    1093            0 :          CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
    1094            0 :          DO gamma = 1, 3
    1095            0 :             DO delta = 1, 3
    1096            0 :                lc_tmp = Levi_Civita(alpha, gamma, delta)
    1097            0 :                IF (lc_tmp == 0._dp) CYCLE
    1098            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
    1099            0 :                CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
    1100            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
    1101            0 :                CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
    1102              : 
    1103              :                ! < mu | nu > * R^mu_gamma
    1104              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
    1105              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
    1106            0 :                                          gauge_origin=gauge_origin)
    1107              : 
    1108              :                ! < mu | nu > * R^nu_gamma
    1109              :                CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
    1110              :                                          qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
    1111            0 :                                          gauge_origin=gauge_origin)
    1112              : 
    1113              :                ! !! A = alpha*A + beta*B or
    1114              :                ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
    1115              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
    1116            0 :                               vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
    1117              : 
    1118              :                ! and accumulate the results
    1119              :                CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
    1120            0 :                               1._dp, vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
    1121              :             END DO
    1122              :          END DO
    1123              : 
    1124              :          CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
    1125            0 :                                       buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
    1126              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
    1127              :                             -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &   ! again, the negative sign is here
    1128            0 :                             1.0_dp, vcd_env%op_dB(ispin))
    1129              : 
    1130              :          ! CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
    1131            0 :          CALL cp_fm_release(buf)
    1132              : 
    1133              :       END ASSOCIATE
    1134              : 
    1135              :       ! We have built op_dB but plug -op_dB into the linres_solver
    1136              :       ! CALL cp_fm_scale(-1.0_dp, vcd_env%op_dB(ispin))
    1137              : 
    1138            0 :       CALL timestop(handle)
    1139            0 :    END SUBROUTINE mfp_build_operator_gauge_independent
    1140              : 
    1141              : ! *****************************************************************************
    1142              : !> \brief Get the dC/dB using the vcd_env%op_dB
    1143              : !> \param vcd_env ...
    1144              : !> \param p_env ...
    1145              : !> \param qs_env ...
    1146              : !> \param alpha ...
    1147              : !> \author Edward Ditler
    1148              : ! **************************************************************************************************
    1149            0 :    SUBROUTINE mfp_response(vcd_env, p_env, qs_env, alpha)
    1150              : 
    1151              :       TYPE(vcd_env_type)                                 :: vcd_env
    1152              :       TYPE(qs_p_env_type)                                :: p_env
    1153              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1154              :       INTEGER, INTENT(IN)                                :: alpha
    1155              : 
    1156              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'mfp_response'
    1157              : 
    1158              :       INTEGER                                            :: handle, output_unit
    1159              :       LOGICAL                                            :: failure, should_stop
    1160            0 :       TYPE(cp_fm_type), DIMENSION(1)                     :: h1_psi0, psi1
    1161              :       TYPE(cp_logger_type), POINTER                      :: logger
    1162              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1163              :       TYPE(linres_control_type), POINTER                 :: linres_control
    1164            0 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1165              :       TYPE(section_vals_type), POINTER                   :: lr_section, vcd_section
    1166              : 
    1167            0 :       CALL timeset(routineN, handle)
    1168            0 :       failure = .FALSE.
    1169              : 
    1170            0 :       NULLIFY (linres_control, lr_section, logger)
    1171              : 
    1172              :       CALL get_qs_env(qs_env=qs_env, &
    1173              :                       dft_control=dft_control, &
    1174              :                       linres_control=linres_control, &
    1175            0 :                       mos=mos)
    1176              : 
    1177            0 :       logger => cp_get_default_logger()
    1178            0 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
    1179              :       vcd_section => section_vals_get_subs_vals(qs_env%input, &
    1180            0 :                                                 "PROPERTIES%LINRES%VCD")
    1181              : 
    1182              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
    1183            0 :                                          extension=".linresLog")
    1184            0 :       IF (output_unit > 0) THEN
    1185              :          WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
    1186            0 :             "*** Self consistent optimization of the magnetic response wavefunction ***"
    1187              :       END IF
    1188              : 
    1189              :       ! allocate the vectors
    1190              :       ASSOCIATE (psi0_order => vcd_env%dcdr_env%mo_coeff)
    1191            0 :          CALL cp_fm_create(psi1(1), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
    1192            0 :          CALL cp_fm_create(h1_psi0(1), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
    1193              : 
    1194              :          ! Restart
    1195            0 :          IF (linres_control%linres_restart) THEN
    1196            0 :             CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, alpha, "dCdB")
    1197              :          ELSE
    1198            0 :             CALL cp_fm_set_all(psi1(1), 0.0_dp)
    1199              :          END IF
    1200              : 
    1201            0 :          IF (output_unit > 0) THEN
    1202              :             WRITE (output_unit, *) &
    1203            0 :                "Response to the perturbation operator referring to the magnetic field in "//ACHAR(alpha + 119)
    1204              :          END IF
    1205              : 
    1206              :          ! First response to get dCR
    1207              :          ! (H0-E0) psi1 = (H1-E1) psi0
    1208              :          ! psi1 = the perturbed wavefunction
    1209              :          ! h1_psi0 = (H1-E1)
    1210              :          ! psi0_order = the unperturbed wavefunction
    1211              :          ! Second response to get dCB
    1212            0 :          CALL cp_fm_set_all(vcd_env%dCB(alpha), 0.0_dp)
    1213            0 :          CALL cp_fm_set_all(h1_psi0(1), 0.0_dp)
    1214            0 :          CALL cp_fm_to_fm(vcd_env%op_dB(1), h1_psi0(1))
    1215              : 
    1216            0 :          linres_control%lr_triplet = .FALSE. ! we do singlet response
    1217            0 :          linres_control%do_kernel = .FALSE. ! no coupled response since imaginary perturbation
    1218            0 :          linres_control%converged = .FALSE.
    1219              :          CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
    1220            0 :                             output_unit, should_stop)
    1221            0 :          CALL cp_fm_to_fm(psi1(1), vcd_env%dCB(alpha))
    1222              : 
    1223              :          ! Write the new result to the restart file
    1224            0 :          IF (linres_control%linres_restart) THEN
    1225            0 :             CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, alpha, "dCdB")
    1226              :          END IF
    1227              : 
    1228              :          ! clean up
    1229            0 :          CALL cp_fm_release(psi1(1))
    1230            0 :          CALL cp_fm_release(h1_psi0(1))
    1231              :       END ASSOCIATE
    1232              : 
    1233              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
    1234            0 :                                         "PRINT%PROGRAM_RUN_INFO")
    1235              : 
    1236            0 :       CALL timestop(handle)
    1237            0 :    END SUBROUTINE mfp_response
    1238              : 
    1239              : ! **************************************************************************************************
    1240              : !> \brief Take matrix < mu | ^O^ | nu > and multiply the blocks with the positions of the basis functions.
    1241              : !>        The matrix consists of blocks (iatom, jatom) corresponding to (mu, nu).
    1242              : !>        With basis_function_nu = .TRUE.  we have to multiply each block by particle_set(jatom)%r(delta)
    1243              : !>        With basis_function_nu = .FALSE. we have to multiply each block by particle_set(iatom)%r(delta)
    1244              : !> \param matrix              The matrix we are operating on
    1245              : !> \param qs_kind_set         Needed to set up the basis_sets
    1246              : !> \param particle_set        Needed for the coordinates
    1247              : !> \param basis_type          Needed to set up the basis_sets
    1248              : !> \param sab_nl              The supplied neighborlist
    1249              : !> \param direction           The index of the nuclear position we are multiplying by
    1250              : !> \param basis_function_nu   Determines whether we multiply by R^nu or R^mu
    1251              : !> \param gauge_origin        The gauge origin, if we do distributed gauge
    1252              : !> \author Edward Ditler
    1253              : ! **************************************************************************************************
    1254            0 :    SUBROUTINE multiply_by_position(matrix, qs_kind_set, particle_set, basis_type, sab_nl, &
    1255              :                                    direction, basis_function_nu, gauge_origin)
    1256              : 
    1257              :       TYPE(dbcsr_type)                                   :: matrix
    1258              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1259              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1260              :       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
    1261              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1262              :          POINTER                                         :: sab_nl
    1263              :       INTEGER                                            :: direction
    1264              :       LOGICAL                                            :: basis_function_nu
    1265              :       REAL(dp), DIMENSION(3), OPTIONAL                   :: gauge_origin
    1266              : 
    1267              :       CHARACTER(len=*), PARAMETER :: routineN = 'multiply_by_position'
    1268              : 
    1269              :       INTEGER                                            :: handle, iatom, icol, ikind, inode, irow, &
    1270              :                                                             jatom, jkind, last_jatom, mepos, &
    1271              :                                                             nkind, nseta, nsetb, nthread
    1272              :       INTEGER, DIMENSION(3)                              :: cell
    1273            0 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
    1274            0 :                                                             npgfb, nsgfa, nsgfb
    1275            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
    1276              :       LOGICAL                                            :: do_symmetric, found
    1277            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
    1278            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: matrix_block, rpgfa, rpgfb, scon_a, &
    1279            0 :                                                             scon_b, zeta, zetb
    1280            0 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
    1281              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
    1282              :       TYPE(neighbor_list_iterator_p_type), &
    1283            0 :          DIMENSION(:), POINTER                           :: nl_iterator
    1284              : 
    1285            0 :       CALL timeset(routineN, handle)
    1286              : 
    1287              :       ! check for symmetry
    1288            0 :       CPASSERT(SIZE(sab_nl) > 0)
    1289            0 :       CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
    1290              : 
    1291              :       ! prepare basis set and coordinates
    1292            0 :       nkind = SIZE(qs_kind_set)
    1293            0 :       ALLOCATE (basis_set_list(nkind))
    1294            0 :       CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
    1295              : 
    1296              :       nthread = 1
    1297            0 : !$    nthread = omp_get_max_threads()
    1298              :       ! Iterate of neighbor list
    1299            0 :       CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
    1300              : 
    1301              : !$OMP PARALLEL DEFAULT(NONE) &
    1302              : !$OMP SHARED (nthread,nl_iterator, do_symmetric) &
    1303              : !$OMP SHARED (matrix,basis_set_list) &
    1304              : !$OMP SHARED (basis_function_nu, direction, particle_set, gauge_origin) &
    1305              : !$OMP PRIVATE (matrix_block,mepos,ikind,jkind,iatom,jatom,cell) &
    1306              : !$OMP PRIVATE (basis_set_a,basis_set_b) &
    1307              : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a) &
    1308              : !$OMP PRIVATE (zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb) &
    1309            0 : !$OMP PRIVATE (zetb, scon_a, scon_b, irow, icol, found, inode, last_jatom)
    1310              :       mepos = 0
    1311              : !$    mepos = omp_get_thread_num()
    1312              : 
    1313              :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
    1314              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
    1315              :                                 iatom=iatom, jatom=jatom, inode=inode)
    1316              :          basis_set_a => basis_set_list(ikind)%gto_basis_set
    1317              :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
    1318              :          basis_set_b => basis_set_list(jkind)%gto_basis_set
    1319              :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
    1320              :          ! basis ikind
    1321              :          ! basis ikind
    1322              :          first_sgfa => basis_set_a%first_sgf
    1323              :          la_max => basis_set_a%lmax
    1324              :          la_min => basis_set_a%lmin
    1325              :          npgfa => basis_set_a%npgf
    1326              :          nseta = basis_set_a%nset
    1327              :          nsgfa => basis_set_a%nsgf_set
    1328              :          rpgfa => basis_set_a%pgf_radius
    1329              :          set_radius_a => basis_set_a%set_radius
    1330              :          scon_a => basis_set_a%scon
    1331              :          zeta => basis_set_a%zet
    1332              :          ! basis jkind
    1333              :          first_sgfb => basis_set_b%first_sgf
    1334              :          lb_max => basis_set_b%lmax
    1335              :          lb_min => basis_set_b%lmin
    1336              :          npgfb => basis_set_b%npgf
    1337              :          nsetb = basis_set_b%nset
    1338              :          nsgfb => basis_set_b%nsgf_set
    1339              :          rpgfb => basis_set_b%pgf_radius
    1340              :          set_radius_b => basis_set_b%set_radius
    1341              :          scon_b => basis_set_b%scon
    1342              :          zetb => basis_set_b%zet
    1343              : 
    1344              :          IF (inode == 1) last_jatom = 0
    1345              : 
    1346              :          ! this guarantees minimum image convention
    1347              :          ! anything else would not make sense
    1348              :          IF (jatom == last_jatom) THEN
    1349              :             CYCLE
    1350              :          END IF
    1351              : 
    1352              :          last_jatom = jatom
    1353              : 
    1354              :          IF (do_symmetric) THEN
    1355              :             IF (iatom <= jatom) THEN
    1356              :                irow = iatom
    1357              :                icol = jatom
    1358              :             ELSE
    1359              :                irow = jatom
    1360              :                icol = iatom
    1361              :             END IF
    1362              :          ELSE
    1363              :             irow = iatom
    1364              :             icol = jatom
    1365              :          END IF
    1366              : 
    1367              :          NULLIFY (matrix_block)
    1368              :          CALL dbcsr_get_block_p(matrix, irow, icol, matrix_block, found)
    1369              :          CPASSERT(found)
    1370              : 
    1371              :          IF (PRESENT(gauge_origin)) THEN
    1372              :             IF (basis_function_nu) THEN
    1373              :                !$OMP CRITICAL(blockadd)
    1374              :                matrix_block(:, :) = matrix_block(:, :)*(particle_set(jatom)%r(direction) - gauge_origin(direction))
    1375              :                !$OMP END CRITICAL(blockadd)
    1376              :             ELSE
    1377              :                !$OMP CRITICAL(blockadd)
    1378              :                matrix_block(:, :) = matrix_block(:, :)*(particle_set(iatom)%r(direction) - gauge_origin(direction))
    1379              :                !$OMP END CRITICAL(blockadd)
    1380              :             END IF
    1381              :          ELSE IF (.NOT. PRESENT(gauge_origin)) THEN
    1382              :             IF (basis_function_nu) THEN
    1383              :                !$OMP CRITICAL(blockadd)
    1384              :                matrix_block(:, :) = matrix_block(:, :)*particle_set(jatom)%r(direction)
    1385              :                !$OMP END CRITICAL(blockadd)
    1386              :             ELSE
    1387              :                !$OMP CRITICAL(blockadd)
    1388              :                matrix_block(:, :) = matrix_block(:, :)*particle_set(iatom)%r(direction)
    1389              :                !$OMP END CRITICAL(blockadd)
    1390              :             END IF
    1391              :          END IF
    1392              :       END DO
    1393              : !$OMP END PARALLEL
    1394            0 :       CALL neighbor_list_iterator_release(nl_iterator)
    1395              : 
    1396              :       ! Release work storage
    1397            0 :       DEALLOCATE (basis_set_list)
    1398              : 
    1399            0 :       CALL timestop(handle)
    1400              : 
    1401            0 :    END SUBROUTINE multiply_by_position
    1402              : 
    1403              : END MODULE qs_mfp
        

Generated by: LCOV version 2.0-1