(git:ed6f26b)
Loading...
Searching...
No Matches
qs_vcd.F
Go to the documentation of this file.
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!--------------------------------------------------------------------------------------------------!
7MODULE qs_vcd
8 USE atomic_kind_types, ONLY: get_atomic_kind
9 USE cell_types, ONLY: cell_type
10 USE commutator_rpnl, ONLY: build_com_mom_nl
11 USE cp_control_types, ONLY: dft_control_type
12 USE cp_dbcsr_api, ONLY: dbcsr_add,&
13 dbcsr_copy,&
14 dbcsr_desymmetrize,&
15 dbcsr_set
16 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
17 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
18 cp_fm_scale_and_add,&
19 cp_fm_trace
20 USE cp_fm_types, ONLY: cp_fm_create,&
21 cp_fm_release,&
22 cp_fm_set_all,&
23 cp_fm_to_fm,&
24 cp_fm_type
25 USE cp_log_handling, ONLY: cp_get_default_logger,&
26 cp_logger_type
27 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
28 cp_print_key_unit_nr
29 USE input_section_types, ONLY: section_vals_get_subs_vals,&
30 section_vals_type
31 USE kinds, ONLY: dp
32 USE parallel_gemm_api, ONLY: parallel_gemm
33 USE particle_types, ONLY: particle_type
34 USE qs_dcdr_ao, ONLY: hr_mult_by_delta_1d
35 USE qs_environment_types, ONLY: get_qs_env,&
36 qs_environment_type
37 USE qs_kind_types, ONLY: get_qs_kind,&
38 qs_kind_type
39 USE qs_linres_methods, ONLY: linres_solver
40 USE qs_linres_types, ONLY: linres_control_type,&
41 vcd_env_type
42 USE qs_mo_types, ONLY: mo_set_type
43 USE qs_moments, ONLY: dipole_velocity_deriv
44 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
45 USE qs_p_env_types, ONLY: qs_p_env_type
46 USE qs_vcd_ao, ONLY: build_dsdv_matrix,&
47 build_dcom_rpnl,&
48 build_drpnl_matrix,&
49 build_matrix_hr_rh,&
50 hr_mult_by_delta_3d
51 USE qs_vcd_utils, ONLY: vcd_read_restart,&
52 vcd_write_restart
53#include "./base/base_uses.f90"
54
55 IMPLICIT NONE
56
57 PRIVATE
58 PUBLIC :: prepare_per_atom_vcd
59 PUBLIC :: vcd_build_op_dv
60 PUBLIC :: vcd_response_dv
61 PUBLIC :: apt_dv
62 PUBLIC :: aat_dv
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd'
65
66 REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = reshape((/ &
67 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
68 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
69 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
70 (/3, 3, 3/))
71 INTEGER, DIMENSION(3, 3), PARAMETER :: multipole_2d_to_1d = reshape([4, 5, 6, 5, 7, 8, 6, 8, 9], [3, 3])
72CONTAINS
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
1126END MODULE qs_vcd
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:75
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:517
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
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:552
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:1665
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:528
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:164
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:853
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:1063
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:731
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, zatom, 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_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_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
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:128
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_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_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_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
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_linres_types::vcd_env_type
Definition qs_linres_types.F:305
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_p_env_types::qs_p_env_type
Represent a qs system that is perturbed. Can calculate the linear operator and the rhs of the system ...
Definition qs_p_env_types.F:58