(git:3add494)
qs_vcd.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 MODULE qs_vcd
8  USE atomic_kind_types, ONLY: get_atomic_kind
9  USE cell_types, ONLY: cell_type
10  USE commutator_rpnl, ONLY: build_com_mom_nl
11  USE cp_control_types, ONLY: dft_control_type
12  USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
13  USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
14  cp_fm_scale_and_add,&
15  cp_fm_trace
16  USE cp_fm_types, ONLY: cp_fm_create,&
17  cp_fm_release,&
18  cp_fm_set_all,&
19  cp_fm_to_fm,&
20  cp_fm_type
21  USE cp_log_handling, ONLY: cp_get_default_logger,&
22  cp_logger_type
23  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
24  cp_print_key_unit_nr
25  USE dbcsr_api, ONLY: dbcsr_add,&
26  dbcsr_copy,&
27  dbcsr_desymmetrize,&
28  dbcsr_set
29  USE input_section_types, ONLY: section_vals_get_subs_vals,&
30  section_vals_type
31  USE kinds, ONLY: dp
32  USE parallel_gemm_api, ONLY: parallel_gemm
33  USE particle_types, ONLY: particle_type
34  USE qs_dcdr_ao, ONLY: hr_mult_by_delta_1d
35  USE qs_environment_types, ONLY: get_qs_env,&
36  qs_environment_type
37  USE qs_kind_types, ONLY: get_qs_kind,&
38  qs_kind_type
39  USE qs_linres_methods, ONLY: linres_solver
40  USE qs_linres_types, ONLY: linres_control_type,&
41  vcd_env_type
42  USE qs_mo_types, ONLY: mo_set_type
43  USE qs_moments, ONLY: dipole_velocity_deriv
44  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
45  USE qs_p_env_types, ONLY: qs_p_env_type
46  USE qs_vcd_ao, ONLY: build_dsdv_matrix,&
47  build_dcom_rpnl,&
48  build_drpnl_matrix,&
49  build_matrix_hr_rh,&
50  hr_mult_by_delta_3d
51  USE qs_vcd_utils, ONLY: vcd_read_restart,&
52  vcd_write_restart
53 #include "./base/base_uses.f90"
54 
55  IMPLICIT NONE
56 
57  PRIVATE
58  PUBLIC :: prepare_per_atom_vcd
59  PUBLIC :: vcd_build_op_dv
60  PUBLIC :: vcd_response_dv
61  PUBLIC :: apt_dv
62  PUBLIC :: aat_dv
63 
64  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd'
65 
66  REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = reshape((/ &
67  0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
68  0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
69  0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
70  (/3, 3, 3/))
71  INTEGER, DIMENSION(3, 3), PARAMETER :: multipole_2d_to_1d = reshape([4, 5, 6, 5, 7, 8, 6, 8, 9], [3, 3])
72 CONTAINS
73 
74 ! **************************************************************************************************
75 !> \brief Compute I_{alpha beta}^lambda = d/dV^lambda_beta <m_alpha> = d/dV^lambda_beta < r x \dot{r} >
76 !> The directions alpha, beta are stored in vcd_env%dcdr_env
77 !> \param vcd_env ...
78 !> \param qs_env ...
79 !> \author Edward Ditler
80 ! **************************************************************************************************
81  SUBROUTINE aat_dv(vcd_env, qs_env)
82  TYPE(vcd_env_type) :: vcd_env
83  TYPE(qs_environment_type), POINTER :: qs_env
84 
85  CHARACTER(LEN=*), PARAMETER :: routinen = 'aat_dV'
86  INTEGER, PARAMETER :: ispin = 1
87 
88  INTEGER :: alpha, delta, gamma, handle, ikind, &
89  my_index, nao, nmo, nspins
90  LOGICAL :: ghost
91  REAL(dp) :: aat_prefactor, aat_tmp, charge, lc_tmp, &
92  tmp_trace
93  REAL(dp), DIMENSION(3, 3) :: aat_tmp_33
94  TYPE(cp_fm_type) :: tmp_aomo
95  TYPE(dft_control_type), POINTER :: dft_control
96  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
97  POINTER :: sab_all, sab_orb, sap_ppnl
98  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
99  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
100 
101  CALL timeset(routinen, handle)
102 
103  CALL get_qs_env(qs_env=qs_env, &
104  dft_control=dft_control, &
105  sap_ppnl=sap_ppnl, &
106  sab_orb=sab_orb, &
107  sab_all=sab_all, &
108  particle_set=particle_set, &
109  qs_kind_set=qs_kind_set)
110 
111  CALL cp_fm_create(tmp_aomo, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
112 
113  nspins = dft_control%nspins
114  nmo = vcd_env%dcdr_env%nmo(ispin)
115  nao = vcd_env%dcdr_env%nao
116  associate(mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), aat_atom => vcd_env%aat_atom_nvpt)
117 
118  ! I_{alpha beta}^lambda = 1/2c \sum_j^occ ...
119  aat_prefactor = 1.0_dp!/(c_light_au * 2._dp)
120  IF (nspins .EQ. 1) aat_prefactor = aat_prefactor*2.0_dp
121 
122  ! The non-PP part of the AAT consists of four contributions:
123  ! (A1): + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda)
124  ! (A2): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta ∂_delta | nu > * (nu == lambda)
125  ! (B): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma | nu > * (delta == beta) * (nu == lambda)
126  ! (C): + iP^1 * ε_{alpha gamma delta} * < mu | r_gamma ∂_delta | nu >
127 
128  ! (A1) + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda)
129  ! (A2) - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta ∂_delta | nu > * (nu == lambda)
130  ! Conjecture : It doesn't matter that the beta and gamma are swapped around!
131  ! We define o = | ∂_delta nu >
132  ! and then < a | r_beta r_gamma | o > = < a | r_gamma r_beta | o>
133  ! (A) + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda - nu == lambda)
134  ! We have built the matrices - < mu | r_beta r_gamma ∂_delta | nu > in vcd_env%moments_der
135  ! moments_der(1:9; 1:3) = moments_der(x, y, z, xx, xy, xz, yy, yz, zz;
136  ! x, y, z)
137 
138  aat_tmp_33 = 0._dp
139  DO gamma = 1, 3
140  my_index = multipole_2d_to_1d(vcd_env%dcdr_env%beta, gamma)
141  DO delta = 1, 3
142  ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
143  ! matrix_nosym_temp = - < mu | r_beta r_gamma ∂_delta | nu > * (mu - nu)
144  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
145  vcd_env%moments_der_right(my_index, delta)%matrix)
146  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
147  vcd_env%moments_der_left(my_index, delta)%matrix, &
148  1._dp, -1._dp)
149 
150  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
151  CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp_33(gamma, delta))
152  END DO
153  END DO
154 
155  DO alpha = 1, 3
156  aat_tmp = 0._dp
157 
158  ! There are two remaining combinations for gamma and delta.
159  DO gamma = 1, 3
160  DO delta = 1, 3
161  lc_tmp = levi_civita(alpha, gamma, delta)
162  IF (lc_tmp == 0._dp) cycle
163 
164  ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
165  ! matrix_nosym_temp = - < mu | r_beta r_gamma ∂_delta | nu > * (mu - nu)
166  ! Because of the negative in moments_der, we need another negative sign here.
167  aat_tmp = aat_tmp + lc_tmp*aat_prefactor*aat_tmp_33(gamma, delta)
168  END DO
169  END DO
170 
171  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
172  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
173  END DO
174 
175  ! (B): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma | nu > * (delta == beta) * (nu == lambda)
176  ! = - P^0 * ε_{alpha gamma beta} * < mu | r_gamma | nu > * (nu == lambda)
177 
178  DO alpha = 1, 3
179  aat_tmp = 0._dp
180 
181  DO gamma = 1, 3
182  lc_tmp = levi_civita(alpha, gamma, vcd_env%dcdr_env%beta)
183  IF (lc_tmp == 0._dp) cycle
184 
185  ! matrix_nosym_temp = < mu | r_gamma | nu > * (nu == lambda)
186  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(gamma)%matrix, &
187  vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
188  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
189  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
190 
191  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
192  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
193  aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace
194  END DO
195 
196  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
197  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
198  END DO
199 
200  ! (C): + iP^1 * ε_{alpha gamma delta} * < mu | r_gamma ∂_delta | nu >
201  DO alpha = 1, 3
202  aat_tmp = 0._dp
203 
204  DO gamma = 1, 3
205  DO delta = 1, 3
206  lc_tmp = levi_civita(alpha, gamma, delta)
207  IF (lc_tmp == 0._dp) cycle
208 
209  CALL cp_dbcsr_sm_fm_multiply(vcd_env%moments_der(gamma, delta)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
210  CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), tmp_trace)
211 
212  ! mo_coeff * dCV_prime = + iP1
213  ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
214  ! so we need the opposite sign.
215  aat_tmp = aat_tmp - 2._dp*aat_prefactor*tmp_trace*lc_tmp
216  END DO
217  END DO
218 
219  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
220  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
221  END DO
222 
223  ! The PP part consists of four contributions
224  ! (D): - P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma [V, r_delta] | nu > * (mu == lambda)
225  ! (E): + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
226  ! (F): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma [[V, r_beta], r_delta] | nu > * (eta == lambda)
227  ! (G): - iP^1 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] | nu >
228 
229  ! (D): - P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma [V, r_delta] | nu > * (mu == lambda)
230  ! The negative of this is in vcd_env%matrix_r_rxvr
231  DO alpha = 1, 3
232  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
233  vcd_env%matrix_r_rxvr(alpha, vcd_env%dcdr_env%beta)%matrix)
234  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
235  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
236 
237  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
238  CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp)
239  aat_tmp = -aat_prefactor*aat_tmp
240 
241  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
242  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
243  END DO
244 
245  ! (E): + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
246  ! This is in vcd_env%matrix_rxvr_r
247  DO alpha = 1, 3
248  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rxvr_r(alpha, vcd_env%dcdr_env%beta)%matrix)
249  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
250  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
251 
252  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
253  CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp)
254  aat_tmp = aat_prefactor*aat_tmp
255 
256  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
257  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
258  END DO
259 
260  ! (F): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma [[V, r_beta], r_delta] | nu > * (eta == lambda)
261  ! + P^0 * ε_{alpha gamma delta} * < mu | [[V, r_beta], r_delta] | nu > * (eta == lambda) * R_gamma
262  ! The negative is in vcd_env%matrix_r_doublecom
263  DO alpha = 1, 3
264  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_r_doublecom(alpha, vcd_env%dcdr_env%beta)%matrix, &
265  mo_coeff, tmp_aomo, ncol=nmo)
266  CALL cp_fm_trace(mo_coeff, tmp_aomo, aat_tmp)
267  aat_tmp = -aat_prefactor*aat_tmp
268 
269  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
270  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
271  END DO
272 
273  ! (G): - iP^1 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] | nu >
274  DO alpha = 1, 3
275  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_rxrv(alpha)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
276  CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), aat_tmp)
277 
278  ! I can take the positive, because build_com_mom_nl computes r x [r, V]
279  aat_tmp = 2._dp*aat_prefactor*aat_tmp
280 
281  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
282  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
283  + aat_tmp
284  END DO
285 
286  ! All the reference dependent stuff
287  ! (C) iP^1 * ε_{alpha gamma delta} * < mu | ∂_delta | nu > * (- R_gamma)
288  DO alpha = 1, 3
289  aat_tmp = 0._dp
290 
291  DO gamma = 1, 3
292  DO delta = 1, 3
293  lc_tmp = levi_civita(alpha, gamma, delta)
294  IF (lc_tmp == 0._dp) cycle
295  ! dipvel_ao = + < a | ∂ | b >
296  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao(delta)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
297  CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), tmp_trace)
298 
299  ! The negative sign is due to (r - O^mag_gamma) and otherwise this is
300  ! exactly the APT dipvel(beta, delta) * (-O^mag_gamma)
301  aat_tmp = aat_tmp + 2._dp*aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
302  END DO
303  END DO
304 
305  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
306  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
307  END DO
308 
309  ! (G): - iP^1 * ε_{alpha gamma delta} * < mu | [V, r_delta] | nu > * (- R_gamma)
310  DO alpha = 1, 3
311  aat_tmp = 0._dp
312  DO gamma = 1, 3
313  DO delta = 1, 3
314  lc_tmp = levi_civita(alpha, gamma, delta)
315  IF (lc_tmp == 0._dp) cycle
316  ! hcom = < a | [r, V] | b > = - < a | [V, r] | b >
317  ! mo_coeff * dCV_prime = + iP1
318  CALL cp_dbcsr_sm_fm_multiply(vcd_env%hcom(delta)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
319  CALL cp_fm_trace(tmp_aomo, vcd_env%dCV_prime(ispin), tmp_trace)
320 
321  ! This is exactly APT hcom(beta, delta)
322  aat_tmp = aat_tmp + 2._dp*aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
323  END DO
324  END DO
325 
326  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
327  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
328  END DO
329 
330  ! mag_vel, vel, mag
331  ! Ai) + ε_{alpha gamma delta} * R_beta R_gamma * < mu | ∂_delta | nu > * (mu - nu)
332  ! Aii) + ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma ∂_delta | nu > * (mu - nu)
333  ! Aiii) + ε_{alpha gamma delta} * (-R_gamma) * < mu | r_beta ∂_delta | nu > * (mu - nu)
334  DO alpha = 1, 3
335  aat_tmp = 0._dp
336  DO gamma = 1, 3
337  DO delta = 1, 3
338  lc_tmp = levi_civita(alpha, gamma, delta)
339  IF (lc_tmp == 0._dp) cycle
340  ! iii) - R_gamma * < mu | r_beta ∂_delta | nu > * (mu - nu)
341  ! mag
342  ! matrix_difdip2(beta, alpha) = - < a | r_beta | ∂_alpha b > * (mu - nu)
343  ! so I need matrix_difdip2(beta, delta)
344  ! Only this part correspond to the APT difdip(beta, alpha)
345  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(vcd_env%dcdr_env%beta, delta)%matrix, mo_coeff, &
346  tmp_aomo, ncol=nmo)
347  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
348 
349  ! There is a negative sign here, because the routine dipole_velocity_deriv calculates
350  ! the derivatives with respect to nuclear positions and we need electronic positions
351  aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace*(-vcd_env%magnetic_origin_atom(gamma))
352 
353  ! This part doesn't appear in the APT
354  ! ii) - R_beta * < mu | r_gamma ∂_delta | nu > * (mu - nu)
355  ! vel
356  ! matrix_difdip2(beta, alpha) = - < a | r_beta | ∂_alpha b > * (mu - nu)
357  ! so I need matrix_difdip2(gamma, delta)
358  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(gamma, delta)%matrix, mo_coeff, &
359  tmp_aomo, ncol=nmo)
360  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
361 
362  ! There is a negative sign here, because the routine dipole_velocity_deriv calculates
363  ! the derivatives with respect to nuclear positions and we need electronic positions
364  aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace*(-vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta))
365 
366  ! i) + R_beta R_gamma * < mu | ∂_delta | nu > * (mu - nu)
367  ! mag_vel
368  ! dipvel_ao = + < a | ∂ | b >
369  CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
370  CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
371  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
372  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
373  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
374  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
375  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
376  1._dp, -1._dp)
377 
378  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
379  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
380  aat_tmp = aat_tmp + lc_tmp*aat_prefactor*tmp_trace* &
381  (vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
382 
383  END DO
384  END DO
385 
386  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
387  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
388  END DO
389 
390  ! (B): P^0 * ε_{alpha gamma beta} * < mu | nu > * (nu == lambda) * R_gamma
391  DO alpha = 1, 3
392  aat_tmp = 0._dp
393 
394  DO gamma = 1, 3
395  lc_tmp = levi_civita(alpha, gamma, vcd_env%dcdr_env%beta)
396  IF (lc_tmp == 0._dp) cycle
397  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, 0.0_dp)
398  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
399  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", sab_all, &
400  vcd_env%dcdr_env%lambda, direction_or=.true.)
401  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
402  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
403 
404  ! This is in total positive because we are calculating
405  ! -1/2c * P * < a | b > * (delta == beta) * (nu == lambda) * (-R_gamma)
406  ! The whole term corresponds to difdip_s
407  aat_tmp = aat_tmp + lc_tmp*aat_prefactor*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
408  END DO
409 
410  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
411  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
412  END DO
413 
414  ! (D): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta [V, r_delta] | nu > * (mu == lambda)
415  ! mag, vel, mag_vel
416  ! Di) - ε_{alpha gamma delta} * (-R_gamma) * < mu | r_beta [V, r_delta] | nu > * (mu == lambda)
417  ! Dii) - ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (mu == lambda)
418  ! Diii) - ε_{alpha gamma delta} * R_beta R_gamma * < mu | [V, r_delta] | nu > * (mu == lambda)
419 
420  DO alpha = 1, 3
421  aat_tmp = 0._dp
422  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
423 
424  DO gamma = 1, 3
425  DO delta = 1, 3
426  lc_tmp = levi_civita(alpha, gamma, delta)
427  IF (lc_tmp == 0._dp) cycle
428  ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
429 
430  ! This corresponds to rcom
431  ! Di) mag
432  ! -(-R_gamma) * < mu | r_beta [V, r_delta] | nu > * (mu == lambda)
433  ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
434  ! so I need vcd_env%matrix_rrcom(delta, beta)
435  ! The multiplication with delta was not done for all directions
436  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
437  vcd_env%matrix_rrcom(delta, vcd_env%dcdr_env%beta)%matrix)
438  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
439  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
440  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
441  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
442  ! The sign is positive in total, because we have the negative coordinate and the whole term was negative
443  aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*vcd_env%magnetic_origin_atom(gamma)
444 
445  ! This doesn't appear in the APT formula
446  ! Dii) vel
447  ! -(-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (mu == lambda)
448  ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
449  ! so I need vcd_env%matrix_rrcom(delta, gamma)
450  ! The multiplication with delta was already done in SUBROUTINE apt_dV
451  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rrcom(delta, gamma)%matrix)
452  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
453  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
454  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
455  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
456  aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)
457 
458  ! Diii) mag_vel
459  ! - R_beta R_gamma * < mu | [V, r_delta] | nu >
460  ! hcom(delta) = - [V, r_delta]
461  CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
462  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
463  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
464  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
465  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
466  ! No need for a negative sign, because hcom already contains the negative sign.
467  aat_tmp = aat_tmp + &
468  aat_prefactor*tmp_trace*lc_tmp &
469  *(vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
470  END DO
471  END DO
472 
473  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
474  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
475  END DO
476 
477  ! (E): + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
478  ! mag, vel, mag_vel
479  ! Ei) + ε_{alpha gamma delta} * (-R_gamma) * < mu | [V, r_delta] r_beta | nu > * (nu == lambda)
480  ! Eii) + ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (nu == lambda)
481  ! Eiii) + ε_{alpha gamma delta} * R_beta R_gamma * < mu | [V, r_delta] | nu > * (nu == lambda)
482  DO alpha = 1, 3
483  aat_tmp = 0._dp
484  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
485 
486  DO gamma = 1, 3
487  DO delta = 1, 3
488  lc_tmp = levi_civita(alpha, gamma, delta)
489  IF (lc_tmp == 0._dp) cycle
490  ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
491  ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
492 
493  ! This corresponds to rcom
494  ! Ei) mag
495  ! (-R_gamma) * < mu | [V, r_delta] r_beta | nu > * (nu == lambda)
496  ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
497  ! so I need vcd_env%matrix_rcomr(delta, beta)
498  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
499  vcd_env%matrix_rcomr(delta, vcd_env%dcdr_env%beta)%matrix)
500  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
501  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
502 
503  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
504  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
505  aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
506 
507  ! This doesn't appear in the APT formula
508  ! E2) vel
509  ! (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (nu == lambda)
510  ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
511  ! so I need vcd_env%matrix_rrcom(delta, gamma)
512  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rrcom(delta, gamma)%matrix)
513  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
514  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
515 
516  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
517  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
518  aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta))
519 
520  ! E3) mag_vel
521  ! R_beta R_gamma * < mu | [V, r_delta] | nu > * (nu == lambda)
522  CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
523  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
524  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
525 
526  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
527  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
528  ! There has to be a minus here, because hcom = [r, V] = - [V, r]
529  aat_tmp = aat_tmp - &
530  aat_prefactor*tmp_trace*lc_tmp* &
531  (vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
532  END DO
533  END DO
534 
535  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
536  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
537  END DO
538 
539  ! (F): - P^0 * ε_{alpha gamma delta} * < mu | [[V, r_beta], r_delta] | nu > * (eta == lambda) * (-R_gamma)
540  ! This corresponds to APT dcom
541  DO alpha = 1, 3
542  aat_tmp = 0._dp
543 
544  DO gamma = 1, 3
545  DO delta = 1, 3
546  lc_tmp = levi_civita(alpha, gamma, delta)
547  IF (lc_tmp == 0._dp) cycle
548  ! vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta) = - < mu | [ [V, r_beta], r_alpha ] | nu >
549  ! so I need matrix_dcom(delta, vcd_env%dcdr_env%beta)
550  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dcom(delta, vcd_env%dcdr_env%beta)%matrix, &
551  mo_coeff, tmp_aomo, ncol=nmo)
552  CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
553  ! matrix_dcom has the negative sign and we include the negative sign of the coordinate
554  aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
555  END DO
556  END DO
557 
558  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
559  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
560  END DO
561 
562  ! Nuclear contribution
563  CALL get_atomic_kind(particle_set(vcd_env%dcdr_env%lambda)%atomic_kind, kind_number=ikind)
564  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
565  IF (.NOT. ghost) THEN
566  DO alpha = 1, 3
567  aat_tmp = 0._dp
568  DO gamma = 1, 3
569  IF (levi_civita(alpha, gamma, vcd_env%dcdr_env%beta) == 0._dp) cycle
570  aat_tmp = aat_tmp + charge &
571  *levi_civita(alpha, gamma, vcd_env%dcdr_env%beta) &
572  *(particle_set(vcd_env%dcdr_env%lambda)%r(gamma) - vcd_env%magnetic_origin_atom(gamma))
573 
574  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
575  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
576  END DO
577  END DO
578  END IF
579  END associate
580 
581  CALL cp_fm_release(tmp_aomo)
582  CALL timestop(handle)
583  END SUBROUTINE aat_dv
584 
585 ! **************************************************************************************************
586 !> \brief Compute E_{alpha beta}^lambda = d/dV^lambda_beta <\mu_alpha> = d/dV^lambda_beta < \dot{r} >
587 !> The directions alpha, beta are stored in vcd_env%dcdr_env
588 !> \param vcd_env ...
589 !> \param qs_env ...
590 !> \author Edward Ditler, Tomas Zimmermann
591 ! **************************************************************************************************
592  SUBROUTINE apt_dv(vcd_env, qs_env)
593  TYPE(vcd_env_type) :: vcd_env
594  TYPE(qs_environment_type), POINTER :: qs_env
595 
596  CHARACTER(LEN=*), PARAMETER :: routinen = 'apt_dV'
597  INTEGER, PARAMETER :: ispin = 1
598  REAL(dp), PARAMETER :: f_spin = 2._dp
599 
600  INTEGER :: alpha, handle, ikind, nao, nmo
601  LOGICAL :: ghost
602  REAL(dp) :: charge
603  REAL(kind=dp) :: apt_dcom, apt_difdip, apt_dipvel, &
604  apt_hcom, apt_rcom
605  TYPE(cp_fm_type) :: buf, matrix_dsdv_mo
606  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
607  POINTER :: sab_all
608  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
609  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
610 
611  CALL timeset(routinen, handle)
612 
613  CALL get_qs_env(qs_env=qs_env, &
614  sab_all=sab_all, &
615  particle_set=particle_set, &
616  qs_kind_set=qs_kind_set)
617 
618  nmo = vcd_env%dcdr_env%nmo(ispin)
619  nao = vcd_env%dcdr_env%nao
620 
621  associate(apt_el => vcd_env%apt_el_nvpt, &
622  apt_nuc => vcd_env%apt_nuc_nvpt, &
623  apt_total => vcd_env%apt_total_nvpt, &
624  mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), &
625  deltar => vcd_env%dcdr_env%deltaR)
626 
627  ! build the full matrices
628  CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
629  CALL cp_fm_create(matrix_dsdv_mo, vcd_env%dcdr_env%momo_fm_struct(ispin)%struct)
630 
631  ! STEP 1: dCV contribution (dipvel + commutator)
632  ! <mu|∂_alpha|nu> and <mu|[r_alpha, V]|nu> in AO basis
633  ! We compute tr(c_1^* x ∂_munu x c_0) + tr(c_0 x ∂_munu x c_1)
634  ! We compute tr(c_1^* x [,]_munu x c_0) + tr(c_0 x [,]_munu x c_1)
635  CALL cp_fm_scale_and_add(0._dp, vcd_env%dCV_prime(ispin), -1._dp, vcd_env%dCV(ispin))
636 
637  ! Ref independent
638  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
639  buf, ncol=nmo)
640  CALL parallel_gemm("T", "N", nmo, nmo, nao, &
641  1.0_dp, mo_coeff, buf, &
642  0.0_dp, matrix_dsdv_mo)
643 
644  CALL parallel_gemm("N", "N", nao, nmo, nmo, &
645  -0.5_dp, mo_coeff, matrix_dsdv_mo, &
646  1.0_dp, vcd_env%dCV_prime(ispin))
647 
648  ! + i∂ - i[Vnl, r]
649  DO alpha = 1, 3
650  CALL cp_fm_set_all(buf, 0.0_dp)
651  apt_dipvel = 0.0_dp
652 
653  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao(alpha)%matrix, mo_coeff, buf, ncol=nmo)
654  CALL cp_fm_trace(buf, vcd_env%dCV_prime(ispin), apt_dipvel)
655  ! dipvel_ao = + < a | ∂ | b >
656  ! mo_coeff * dCV_prime = + iP1
657  apt_dipvel = 2._dp*apt_dipvel
658  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
659  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_dipvel
660  END DO
661 
662  DO alpha = 1, 3
663  CALL cp_fm_set_all(buf, 0.0_dp)
664  apt_hcom = 0.0_dp
665  CALL cp_dbcsr_sm_fm_multiply(vcd_env%hcom(alpha)%matrix, mo_coeff, buf, ncol=nmo)
666  CALL cp_fm_trace(buf, vcd_env%dCV_prime(ispin), apt_hcom)
667 
668  ! hcom = < a | [r, V] | b > = - < a | [V, r] | b >
669  ! mo_coeff * dCV_prime = + iP1
670  apt_hcom = +2._dp*apt_hcom
671 
672  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
673  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_hcom
674  END DO !x/y/z
675 
676  ! STEP 2: basis function derivative contribution
677  !! difdip_s
678  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, 0.0_dp)
679  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s1(1)%matrix, &
680  vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix)
681  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, qs_kind_set, "ORB", sab_all, &
682  vcd_env%dcdr_env%lambda, direction_or=.true.)
683 
684  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
685  buf, ncol=nmo, alpha=1._dp, beta=0._dp)
686  CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
687 
688  apt_difdip = -f_spin*apt_difdip
689  apt_el(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) &
690  = apt_el(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) + apt_difdip
691 
692  !! difdip(j, idir) = < a | r_j | ∂_idir b >
693  !! matrix_difdip2(beta, alpha) = < a | r_beta | ∂_alpha b >
694  DO alpha = 1, 3 ! x/y/z for differentiated AO
695  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(vcd_env%dcdr_env%beta, alpha)%matrix, mo_coeff, &
696  buf, ncol=nmo, alpha=1._dp, beta=0._dp)
697 
698  CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
699  ! There is a negative sign here, because the routine dipole_velocity_deriv calculates
700  ! the derivatives with respect to nuclear positions and we need electronic positions
701  apt_difdip = -f_spin*apt_difdip
702  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
703  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + apt_difdip
704 
705  END DO !alpha
706 
707  ! STEP 3: The terms r * [V, r]
708  ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
709  ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
710  DO alpha = 1, 3 ! x/y/z for differentiated AO
711  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rcomr(alpha, vcd_env%dcdr_env%beta)%matrix)
712  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rrcom(alpha, vcd_env%dcdr_env%beta)%matrix)
713 
714  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
715  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
716  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
717  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
718 
719  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
720  1.0_dp, -1.0_dp)
721 
722  CALL cp_fm_set_all(buf, 0.0_dp)
723  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
724  CALL cp_fm_trace(mo_coeff, buf, apt_rcom)
725 
726  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
727  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_rcom
728  END DO !alpha
729 
730  ! STEP 4: pseudopotential derivative contribution
731  ! vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta) = - < mu | [ [V, r_beta], r_alpha ] | nu >
732  DO alpha = 1, 3 !x/y/z for differentiated AO
733  CALL cp_fm_set_all(buf, 0.0_dp)
734  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta)%matrix, mo_coeff, buf, ncol=nmo)
735  CALL cp_fm_trace(mo_coeff, buf, apt_dcom)
736  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
737  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_dcom
738  END DO !alpha
739 
740  ! The reference point dependent terms:
741  !! difdip_munu
742  ! The additional term here is < a | db/dr(alpha)> * (delta_a - delta_b) * ref_point(beta)
743  ! in qs_env%matrix_s1(2:4) there is < da/dR | b > = - < da/dr | b > = < a | db/dr >
744  DO alpha = 1, 3
745  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, 0._dp)
746  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, 0._dp)
747  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(alpha + 1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
748  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
749 
750  ! < a | db/dr(alpha) > * R^lambda_beta * delta^lambda_nu
751  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
752  vcd_env%dcdr_env%lambda, direction_or=.true.)
753  ! < a | db/dr(alpha) > * R^lambda_beta * delta^lambda_mu
754  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
755  vcd_env%dcdr_env%lambda, direction_or=.false.)
756 
757  ! < a | db/dr > * R^lambda_beta * ( delta^lambda_mu - delta^lambda_nu )
758  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, &
759  1._dp, -1._dp)
760 
761  CALL cp_fm_set_all(buf, 0.0_dp)
762  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, mo_coeff, buf, ncol=nmo)
763  CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
764 
765  ! And the whole contribution is
766  ! - < a | db/dr > * (mu - nu) * ref_point
767  apt_difdip = -apt_difdip*vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)
768 
769  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
770  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_difdip
771  END DO
772 
773  ! And the additional factor to rcom
774  ! < mu | [V, r] | nu > * R^lambda_beta * delta^lambda_mu
775  ! - < mu | [V, r] | nu > * R^lambda_beta * delta^lambda_nu
776  !
777  ! vcd_env%hcom(alpha) = - < mu | [V, r_alpha] | nu >
778  ! particle_set(lambda)%r(vcd_env%dcdr_env%beta) = R^lambda_beta
779 
780  DO alpha = 1, 3
781  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, 0._dp)
782  CALL dbcsr_desymmetrize(vcd_env%hcom(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
783  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
784 
785  ! < mu | [V, r] | nu > * delta^lambda_nu
786  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
787  vcd_env%dcdr_env%lambda, direction_or=.true.)
788  ! < mu | [V, r] | nu > * delta^lambda_mu
789  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
790  vcd_env%dcdr_env%lambda, direction_or=.false.)
791 
792  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, &
793  vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, -1._dp, +1._dp)
794 
795  CALL cp_fm_set_all(buf, 0.0_dp)
796  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, mo_coeff, buf, ncol=nmo)
797  CALL cp_fm_trace(mo_coeff, buf, apt_rcom)
798  apt_rcom = -vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*apt_rcom
799 
800  apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
801  = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_rcom
802  END DO
803 
804  ! STEP 5: nuclear contribution
805  associate(atomic_kind => particle_set(vcd_env%dcdr_env%lambda)%atomic_kind)
806  CALL get_atomic_kind(atomic_kind, kind_number=ikind)
807  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
808  IF (.NOT. ghost) THEN
809  apt_nuc(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) = &
810  apt_nuc(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) + charge
811  END IF
812  END associate
813 
814  ! STEP 6: deallocations
815  CALL cp_fm_release(buf)
816  CALL cp_fm_release(matrix_dsdv_mo)
817 
818  END associate
819 
820  CALL timestop(handle)
821  END SUBROUTINE apt_dv
822 
823 ! **************************************************************************************************
824 !> \brief Initialize the matrices for the NVPT calculation
825 !> \param vcd_env ...
826 !> \param qs_env ...
827 !> \author Edward Ditler
828 ! **************************************************************************************************
829  SUBROUTINE prepare_per_atom_vcd(vcd_env, qs_env)
830  TYPE(vcd_env_type) :: vcd_env
831  TYPE(qs_environment_type), POINTER :: qs_env
832 
833  CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_per_atom_vcd'
834 
835  INTEGER :: handle, i, ispin, j
836  TYPE(cell_type), POINTER :: cell
837  TYPE(dft_control_type), POINTER :: dft_control
838  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
839  POINTER :: sab_all, sab_orb, sap_ppnl
840  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
841  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
842 
843  CALL timeset(routinen, handle)
844 
845  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
846  sab_orb=sab_orb, sab_all=sab_all, sap_ppnl=sap_ppnl, &
847  qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
848 
849  IF (vcd_env%distributed_origin) THEN
850  vcd_env%magnetic_origin_atom(:) = particle_set(vcd_env%dcdr_env%lambda)%r(:) - vcd_env%magnetic_origin(:)
851  vcd_env%spatial_origin_atom = particle_set(vcd_env%dcdr_env%lambda)%r(:) - vcd_env%spatial_origin(:)
852  END IF
853 
854  ! Reset the matrices
855  DO ispin = 1, dft_control%nspins
856  DO j = 1, 3
857  CALL dbcsr_set(vcd_env%matrix_dSdV(j)%matrix, 0._dp)
858  CALL dbcsr_set(vcd_env%matrix_drpnl(j)%matrix, 0._dp)
859 
860  DO i = 1, 3
861  CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0.0_dp)
862  CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0._dp)
863  END DO
864  END DO
865  CALL cp_fm_set_all(vcd_env%op_dV(ispin), 0._dp)
866  CALL dbcsr_set(vcd_env%matrix_hxc_dsdv(ispin)%matrix, 0._dp)
867  END DO
868 
869  ! operator dV
870  ! <mu|d/dV_beta [V, r_alpha]|nu>
871  CALL build_dcom_rpnl(vcd_env%matrix_dcom, qs_kind_set, sab_orb, sap_ppnl, &
872  dft_control%qs_control%eps_ppnl, particle_set, vcd_env%dcdr_env%lambda)
873 
874  ! PP derivative
875  CALL build_drpnl_matrix(vcd_env%matrix_drpnl, qs_kind_set, sab_all, sap_ppnl, &
876  dft_control%qs_control%eps_ppnl, particle_set, pseudoatom=vcd_env%dcdr_env%lambda)
877  ! lin_mom
878  DO i = 1, 3
879  CALL dbcsr_set(vcd_env%dipvel_ao_delta(i)%matrix, 0._dp)
880  CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dipvel_ao(i)%matrix)
881  END DO
882 
883  CALL hr_mult_by_delta_3d(vcd_env%dipvel_ao_delta, qs_kind_set, "ORB", &
884  sab_all, vcd_env%dcdr_env%delta_basis_function, direction_or=.true.)
885 
886  ! dS/dV
887  CALL build_dsdv_matrix(qs_env, vcd_env%matrix_dSdV, &
888  deltar=vcd_env%dcdr_env%delta_basis_function, &
889  rcc=vcd_env%spatial_origin_atom)
890 
891  CALL dipole_velocity_deriv(qs_env, vcd_env%matrix_difdip2, 1, lambda=vcd_env%dcdr_env%lambda, &
892  rc=[0._dp, 0._dp, 0._dp])
893  ! AAT
894  ! moments_throw: x, y, z, xx, xy, xz, yy, yz, zz
895  ! moments_der: (moment, xyz derivative)
896  ! build_local_moments_der_matrix uses adbdr for calculating derivatives of the *primitive*
897  ! on the right. So the resulting
898  ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
899  DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
900  DO j = 1, 3
901  CALL dbcsr_set(vcd_env%moments_der_right(i, j)%matrix, 0.0_dp)
902  CALL dbcsr_set(vcd_env%moments_der_left(i, j)%matrix, 0.0_dp)
903  END DO
904  END DO
905 
906  DO i = 1, 9
907  DO j = 1, 3 ! derivatives
908  CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_right(i, j)%matrix) ! A2
909  CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_left(i, j)%matrix) ! A1
910 
911  ! - < mu | r_beta r_gamma ∂_delta | nu > * (mu/nu == lambda)
912  CALL hr_mult_by_delta_1d(vcd_env%moments_der_right(i, j)%matrix, qs_kind_set, "ORB", &
913  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
914  CALL hr_mult_by_delta_1d(vcd_env%moments_der_left(i, j)%matrix, qs_kind_set, "ORB", &
915  sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
916  END DO
917  END DO
918 
919  DO i = 1, 3
920  DO j = 1, 3
921  CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
922  END DO
923  END DO
924 
925  CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
926  particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
927  matrix_r_doublecom=vcd_env%matrix_r_doublecom, &
928  pseudoatom=vcd_env%dcdr_env%lambda)
929 
930  CALL timestop(handle)
931 
932  END SUBROUTINE prepare_per_atom_vcd
933 
934 ! **************************************************************************************************
935 !> \brief What we are building here is the operator for the NVPT response:
936 !> H0 * C1 - S0 * E0 * C1 = - op_dV
937 !> linres_solver = - [ H1 * C0 - S1 * C0 * E0 ]
938 !> with
939 !> H1 * C0 = dH/dV * C0
940 !> + i[∂]δ * C0
941 !> - i S0 * C^(1,R)
942 !> + i S0 * C0 * (C0 * S^(1,R) * C0)
943 !> - S1 * C0 * E0
944 !>
945 !> H1 * C0 = + i (Hr - rH) * C0 [STEP 1]
946 !> + i[∂]δ * C0 [STEP 2]
947 !> - i[V, r]δ * C0 [STEP 3]
948 !> - i S0 * C^(1,R) [STEP 4]
949 !> - S1 * C0 * E0 [STEP 5]
950 !> \param vcd_env ...
951 !> \param qs_env ...
952 !> \author Edward Ditler, Tomas Zimmermann
953 ! **************************************************************************************************
954  SUBROUTINE vcd_build_op_dv(vcd_env, qs_env)
955  TYPE(vcd_env_type) :: vcd_env
956  TYPE(qs_environment_type), POINTER :: qs_env
957 
958  CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_build_op_dV'
959  INTEGER, PARAMETER :: ispin = 1
960 
961  INTEGER :: handle, nao, nmo
962  TYPE(cp_fm_type) :: buf
963  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
964  POINTER :: sab_all
965  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
966 
967  CALL timeset(routinen, handle)
968 
969  CALL get_qs_env(qs_env=qs_env, &
970  sab_all=sab_all, &
971  qs_kind_set=qs_kind_set)
972 
973  nmo = vcd_env%dcdr_env%nmo(1)
974  nao = vcd_env%dcdr_env%nao
975 
976  CALL build_matrix_hr_rh(vcd_env, qs_env, vcd_env%spatial_origin_atom)
977 
978  ! STEP 1: hr-rh
979  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_hr(ispin, vcd_env%dcdr_env%beta)%matrix)
980  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rh(ispin, vcd_env%dcdr_env%beta)%matrix)
981 
982  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
983  sab_all, vcd_env%dcdr_env%lambda, direction_or=.true.)
984  CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
985  sab_all, vcd_env%dcdr_env%lambda, direction_or=.false.)
986  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
987  vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
988  1.0_dp, -1.0_dp)
989 
990  associate(mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
991  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, &
992  vcd_env%op_dV(ispin), ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
993 
994  ! STEP 2: electronic momentum operator contribution
995  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao_delta(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
996  vcd_env%op_dV(ispin), &
997  ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
998 
999  ! STEP 3: +dV_ppnl/dV, but build_drpnl_matrix gives the negative of dV_ppnl
1000  ! The arguments (-1, 1) are swapped wrt to the hr-rh term, implying that
1001  ! direction_Or and direction_hr do what they should.
1002  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_drpnl(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
1003  vcd_env%op_dV(ispin), &
1004  ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
1005 
1006  ! STEP 4: - S0 * C^(1,R)
1007  CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1008  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), &
1009  vcd_env%op_dV(1), ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
1010 
1011  ! STEP 5: -S(1,V) * C0 * E0
1012  CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
1013  buf, nmo, alpha=1.0_dp, beta=0.0_dp)
1014  CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
1015  -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &
1016  1.0_dp, vcd_env%op_dV(ispin))
1017 
1018  CALL cp_fm_release(buf)
1019  END associate
1020 
1021  ! We have built op_dV but plug -op_dV into the linres_solver
1022  CALL cp_fm_scale(-1.0_dp, vcd_env%op_dV(1))
1023 
1024  ! Revert the matrices
1025  CALL build_matrix_hr_rh(vcd_env, qs_env, [0._dp, 0._dp, 0._dp])
1026 
1027  CALL timestop(handle)
1028  END SUBROUTINE vcd_build_op_dv
1029 
1030 ! *****************************************************************************
1031 !> \brief Get the dC/dV using the vcd_env%op_dV
1032 !> \param vcd_env ...
1033 !> \param p_env ...
1034 !> \param qs_env ...
1035 !> \author Edward Ditler, Tomas Zimmermann
1036 ! **************************************************************************************************
1037  SUBROUTINE vcd_response_dv(vcd_env, p_env, qs_env)
1038 
1039  TYPE(vcd_env_type) :: vcd_env
1040  TYPE(qs_p_env_type) :: p_env
1041  TYPE(qs_environment_type), POINTER :: qs_env
1042 
1043  CHARACTER(LEN=*), PARAMETER :: routinen = 'vcd_response_dV'
1044  INTEGER, PARAMETER :: ispin = 1
1045 
1046  INTEGER :: handle, output_unit
1047  LOGICAL :: failure, should_stop
1048  TYPE(cp_fm_type), DIMENSION(1) :: h1_psi0, psi1
1049  TYPE(cp_logger_type), POINTER :: logger
1050  TYPE(linres_control_type), POINTER :: linres_control
1051  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1052  TYPE(section_vals_type), POINTER :: lr_section, vcd_section
1053 
1054  CALL timeset(routinen, handle)
1055  failure = .false.
1056 
1057  NULLIFY (linres_control, lr_section, logger)
1058 
1059  CALL get_qs_env(qs_env=qs_env, &
1060  linres_control=linres_control, &
1061  mos=mos)
1062 
1063  logger => cp_get_default_logger()
1064  lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
1065  vcd_section => section_vals_get_subs_vals(qs_env%input, &
1066  "PROPERTIES%LINRES%vcd")
1067 
1068  output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
1069  extension=".linresLog")
1070  IF (output_unit > 0) THEN
1071  WRITE (unit=output_unit, fmt="(T10,A,/)") &
1072  "*** Self consistent optimization of the response wavefunction ***"
1073  END IF
1074 
1075  associate(psi0_order => vcd_env%dcdr_env%mo_coeff)
1076  CALL cp_fm_create(psi1(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1077  CALL cp_fm_create(h1_psi0(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1078 
1079  ! Restart
1080  IF (linres_control%linres_restart) THEN
1081  CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
1082  ELSE
1083  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
1084  END IF
1085 
1086  IF (output_unit > 0) THEN
1087  WRITE (output_unit, *) &
1088  "Response to the perturbation operator referring to the velocity of atom ", &
1089  vcd_env%dcdr_env%lambda, " in "//achar(vcd_env%dcdr_env%beta + 119)
1090  END IF
1091 
1092  ! First response to get dCR
1093  ! (H0-E0) psi1 = (H1-E1) psi0
1094  ! psi1 = the perturbed wavefunction
1095  ! h1_psi0 = (H1-E1)
1096  ! psi0_order = the unperturbed wavefunction
1097  ! Second response to get dCV
1098  CALL cp_fm_set_all(vcd_env%dCV(ispin), 0.0_dp)
1099  CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
1100  CALL cp_fm_to_fm(vcd_env%op_dV(ispin), h1_psi0(ispin))
1101 
1102  linres_control%lr_triplet = .false. ! we do singlet response
1103  linres_control%do_kernel = .false. ! no coupled response since imaginary perturbation
1104  linres_control%converged = .false.
1105  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
1106  output_unit, should_stop)
1107  CALL cp_fm_to_fm(psi1(ispin), vcd_env%dCV(ispin))
1108 
1109  ! Write the new result to the restart file
1110  IF (linres_control%linres_restart) THEN
1111  CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
1112  END IF
1113 
1114  END associate
1115 
1116  ! clean up
1117  CALL cp_fm_release(psi1(ispin))
1118  CALL cp_fm_release(h1_psi0(ispin))
1119 
1120  CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1121  "PRINT%PROGRAM_RUN_INFO")
1122 
1123  CALL timestop(handle)
1124  END SUBROUTINE vcd_response_dv
1125 
1126 END MODULE qs_vcd
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition: commutator_rpnl.F:439
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition: cp_fm_basic_linalg.F:1680
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_dcdr_ao
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition: qs_dcdr_ao.F:13
qs_dcdr_ao::hr_mult_by_delta_1d
subroutine, public hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_Or)
Enforce that one of the basis functions in < a | O | b > is centered on atom lambda.
Definition: qs_dcdr_ao.F:703
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_linres_methods
localize wavefunctions linear response scf
Definition: qs_linres_methods.F:15
qs_linres_methods::linres_solver
subroutine, public linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
scf loop to optimize the first order wavefunctions (psi1) given a perturbation as an operator applied...
Definition: qs_linres_methods.F:216
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
qs_moments::dipole_velocity_deriv
subroutine, public dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the basis function on ...
Definition: qs_moments.F:126
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition: qs_p_env_types.F:14
qs_vcd_ao
Definition: qs_vcd_ao.F:7
qs_vcd_ao::build_dcom_rpnl
subroutine, public build_dcom_rpnl(matrix_rv, qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, pseudoatom)
Calculate the double commutator [[Vnl, r], r].
Definition: qs_vcd_ao.F:1854
qs_vcd_ao::hr_mult_by_delta_3d
subroutine, public hr_mult_by_delta_3d(matrix_hr, qs_kind_set, basis_type, sab_nl, deltaR, direction_Or)
Apply the operator \delta_\mu^\lambda to zero out all elements of the matrix which don't fulfill the ...
Definition: qs_vcd_ao.F:2719
qs_vcd_ao::build_matrix_hr_rh
subroutine, public build_matrix_hr_rh(vcd_env, qs_env, rc)
Build the matrix Hr*delta_nu^\lambda - rH*delta_mu^\lambda.
Definition: qs_vcd_ao.F:118
qs_vcd_ao::build_dsdv_matrix
subroutine, public build_dsdv_matrix(qs_env, matrix_dsdv, deltaR, rcc)
Builds the overlap derivative wrt nuclear velocities dS/dV = < mu | r | nu > * (nu - mu)
Definition: qs_vcd_ao.F:1467
qs_vcd_ao::build_drpnl_matrix
subroutine, public build_drpnl_matrix(matrix_rv, qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, pseudoatom)
dV_nl/dV. Adapted from build_com_rpnl.
Definition: qs_vcd_ao.F:2270
qs_vcd_utils
Definition: qs_vcd_utils.F:8
qs_vcd_utils::vcd_write_restart
subroutine, public vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_write_restart.
Definition: qs_vcd_utils.F:658
qs_vcd_utils::vcd_read_restart
subroutine, public vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_read_restart.
Definition: qs_vcd_utils.F:508
qs_vcd
Definition: qs_vcd.F:7
qs_vcd::vcd_build_op_dv
subroutine, public vcd_build_op_dv(vcd_env, qs_env)
What we are building here is the operator for the NVPT response: H0 * C1 - S0 * E0 * C1 = - op_dV lin...
Definition: qs_vcd.F:955
qs_vcd::vcd_response_dv
subroutine, public vcd_response_dv(vcd_env, p_env, qs_env)
Get the dC/dV using the vcd_envop_dV.
Definition: qs_vcd.F:1038
qs_vcd::aat_dv
subroutine, public aat_dv(vcd_env, qs_env)
Compute I_{alpha beta}^lambda = d/dV^lambda_beta <m_alpha> = d/dV^lambda_beta < r x.
Definition: qs_vcd.F:82
qs_vcd::apt_dv
subroutine, public apt_dv(vcd_env, qs_env)
Compute E_{alpha beta}^lambda = d/dV^lambda_beta <\mu_alpha> = d/dV^lambda_beta <.
Definition: qs_vcd.F:593
qs_vcd::prepare_per_atom_vcd
subroutine, public prepare_per_atom_vcd(vcd_env, qs_env)
Initialize the matrices for the NVPT calculation.
Definition: qs_vcd.F:830