(git:f2099e5)
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-2026 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_scale,&
16 dbcsr_set
17 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
18 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
19 cp_fm_scale_and_add,&
20 cp_fm_trace
21 USE cp_fm_types, ONLY: cp_fm_create,&
22 cp_fm_release,&
23 cp_fm_set_all,&
24 cp_fm_to_fm,&
25 cp_fm_type
26 USE cp_log_handling, ONLY: cp_get_default_logger,&
27 cp_logger_type
28 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 cp_print_key_unit_nr
30 USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 section_vals_type
32 USE kinds, ONLY: dp
33 USE parallel_gemm_api, ONLY: parallel_gemm
34 USE particle_types, ONLY: particle_type
35 USE qs_dcdr_ao, ONLY: hr_mult_by_delta_1d
36 USE qs_environment_types, ONLY: get_qs_env,&
37 qs_environment_type
38 USE qs_kind_types, ONLY: get_qs_kind,&
39 qs_kind_type
40 USE qs_linres_methods, ONLY: linres_solver
41 USE qs_linres_types, ONLY: linres_control_type,&
42 vcd_env_type
43 USE qs_mo_types, ONLY: mo_set_type
44 USE qs_moments, ONLY: build_local_moments_der_matrix
45 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
46 USE qs_p_env_types, ONLY: qs_p_env_type
47 USE qs_vcd_ao, ONLY: build_dsdv_matrix,&
48 build_dcom_rpnl,&
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 == 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 ! matrix_difdip2 stores nuclear derivatives; the electronic-coordinate
332 ! derivative contributions below therefore use a negative sign.
333 ! Ai) + ε_{alpha gamma delta} * R_beta R_gamma * < mu | ∂_delta | nu > * (mu - nu)
334 ! Aii) + ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma ∂_delta | nu > * (mu - nu)
335 ! Aiii) + ε_{alpha gamma delta} * (-R_gamma) * < mu | r_beta ∂_delta | nu > * (mu - nu)
336 DO alpha = 1, 3
337 aat_tmp = 0._dp
338 DO gamma = 1, 3
339 DO delta = 1, 3
340 lc_tmp = levi_civita(alpha, gamma, delta)
341 IF (lc_tmp == 0._dp) cycle
342 ! iii) - R_gamma * < mu | r_beta ∂_delta | nu > * (mu - nu)
343 ! mag
344 ! matrix_difdip2(beta, alpha) = - < a | r_beta | ∂_alpha b > * (mu - nu)
345 ! so I need matrix_difdip2(beta, delta)
346 ! Only this part correspond to the APT difdip(beta, alpha)
347 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(vcd_env%dcdr_env%beta, delta)%matrix, mo_coeff, &
348 tmp_aomo, ncol=nmo)
349 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
350
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 aat_tmp = aat_tmp - lc_tmp*aat_prefactor*tmp_trace*(-vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta))
363
364 ! i) + R_beta R_gamma * < mu | ∂_delta | nu > * (mu - nu)
365 ! mag_vel
366 ! dipvel_ao = + < a | ∂ | b >
367 CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
368 CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
369 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
370 sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
371 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
372 sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
373 CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
374 1._dp, -1._dp)
375
376 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
377 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
378 aat_tmp = aat_tmp + lc_tmp*aat_prefactor*tmp_trace* &
379 (vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
380
381 END DO
382 END DO
383
384 aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
385 = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
386 END DO
387
388 ! (B): P^0 * ε_{alpha gamma beta} * < mu | nu > * (nu == lambda) * R_gamma
389 DO alpha = 1, 3
390 aat_tmp = 0._dp
391
392 DO gamma = 1, 3
393 lc_tmp = levi_civita(alpha, gamma, vcd_env%dcdr_env%beta)
394 IF (lc_tmp == 0._dp) cycle
395 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, 0.0_dp)
396 CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
397 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", sab_all, &
398 vcd_env%dcdr_env%lambda, direction_or=.true.)
399 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
400 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
401
402 ! This is in total positive because we are calculating
403 ! -1/2c * P * < a | b > * (delta == beta) * (nu == lambda) * (-R_gamma)
404 ! The whole term corresponds to difdip_s
405 aat_tmp = aat_tmp + lc_tmp*aat_prefactor*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
406 END DO
407
408 aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
409 = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
410 END DO
411
412 ! (D): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta [V, r_delta] | nu > * (mu == lambda)
413 ! mag, vel, mag_vel
414 ! Di) - ε_{alpha gamma delta} * (-R_gamma) * < mu | r_beta [V, r_delta] | nu > * (mu == lambda)
415 ! Dii) - ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (mu == lambda)
416 ! Diii) - ε_{alpha gamma delta} * R_beta R_gamma * < mu | [V, r_delta] | nu > * (mu == lambda)
417
418 DO alpha = 1, 3
419 aat_tmp = 0._dp
420 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
421
422 DO gamma = 1, 3
423 DO delta = 1, 3
424 lc_tmp = levi_civita(alpha, gamma, delta)
425 IF (lc_tmp == 0._dp) cycle
426 ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
427
428 ! This corresponds to rcom
429 ! Di) mag
430 ! -(-R_gamma) * < mu | r_beta [V, r_delta] | nu > * (mu == lambda)
431 ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
432 ! so I need vcd_env%matrix_rrcom(delta, beta)
433 ! The multiplication with delta was not done for all directions
434 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
435 vcd_env%matrix_rrcom(delta, vcd_env%dcdr_env%beta)%matrix)
436 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
437 sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
438 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
439 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
440 ! The sign is positive in total, because we have the negative coordinate and the whole term was negative
441 aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*vcd_env%magnetic_origin_atom(gamma)
442
443 ! This doesn't appear in the APT formula
444 ! Dii) vel
445 ! -(-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (mu == lambda)
446 ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
447 ! so I need vcd_env%matrix_rrcom(delta, gamma)
448 ! The multiplication with delta was already done in SUBROUTINE apt_dV
449 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rrcom(delta, gamma)%matrix)
450 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
451 sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
452 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
453 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
454 aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)
455
456 ! Diii) mag_vel
457 ! - R_beta R_gamma * < mu | [V, r_delta] | nu >
458 ! hcom(delta) = - [V, r_delta]
459 CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
460 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
461 sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
462 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
463 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
464 ! No need for a negative sign, because hcom already contains the negative sign.
465 aat_tmp = aat_tmp + &
466 aat_prefactor*tmp_trace*lc_tmp &
467 *(vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
468 END DO
469 END DO
470
471 aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
472 = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
473 END DO
474
475 ! (E): + P^0 * ε_{alpha gamma delta} * < mu | r_gamma [V, r_delta] r_beta | nu > * (nu == lambda)
476 ! mag, vel, mag_vel
477 ! Ei) + ε_{alpha gamma delta} * (-R_gamma) * < mu | [V, r_delta] r_beta | nu > * (nu == lambda)
478 ! Eii) + ε_{alpha gamma delta} * (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (nu == lambda)
479 ! Eiii) + ε_{alpha gamma delta} * R_beta R_gamma * < mu | [V, r_delta] | nu > * (nu == lambda)
480 DO alpha = 1, 3
481 aat_tmp = 0._dp
482 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
483
484 DO gamma = 1, 3
485 DO delta = 1, 3
486 lc_tmp = levi_civita(alpha, gamma, delta)
487 IF (lc_tmp == 0._dp) cycle
488 ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
489 ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
490
491 ! This corresponds to rcom
492 ! Ei) mag
493 ! (-R_gamma) * < mu | [V, r_delta] r_beta | nu > * (nu == lambda)
494 ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
495 ! so I need vcd_env%matrix_rcomr(delta, beta)
496 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
497 vcd_env%matrix_rcomr(delta, vcd_env%dcdr_env%beta)%matrix)
498 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
499 sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
500
501 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
502 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
503 aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
504
505 ! This doesn't appear in the APT formula
506 ! E2) vel
507 ! (-R_beta) * < mu | r_gamma [V, r_delta] | nu > * (nu == lambda)
508 ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
509 ! so I need vcd_env%matrix_rrcom(delta, gamma)
510 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rrcom(delta, gamma)%matrix)
511 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
512 sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
513
514 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
515 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
516 aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta))
517
518 ! E3) mag_vel
519 ! R_beta R_gamma * < mu | [V, r_delta] | nu > * (nu == lambda)
520 CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
521 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
522 sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
523
524 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, tmp_aomo, ncol=nmo)
525 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
526 ! There has to be a minus here, because hcom = [r, V] = - [V, r]
527 aat_tmp = aat_tmp - &
528 aat_prefactor*tmp_trace*lc_tmp* &
529 (vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*vcd_env%magnetic_origin_atom(gamma))
530 END DO
531 END DO
532
533 aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
534 = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
535 END DO
536
537 ! (F): - P^0 * ε_{alpha gamma delta} * < mu | [[V, r_beta], r_delta] | nu > * (eta == lambda) * (-R_gamma)
538 ! This corresponds to APT dcom
539 DO alpha = 1, 3
540 aat_tmp = 0._dp
541
542 DO gamma = 1, 3
543 DO delta = 1, 3
544 lc_tmp = levi_civita(alpha, gamma, delta)
545 IF (lc_tmp == 0._dp) cycle
546 ! vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta) = - < mu | [ [V, r_beta], r_alpha ] | nu >
547 ! so I need matrix_dcom(delta, vcd_env%dcdr_env%beta)
548 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dcom(delta, vcd_env%dcdr_env%beta)%matrix, &
549 mo_coeff, tmp_aomo, ncol=nmo)
550 CALL cp_fm_trace(mo_coeff, tmp_aomo, tmp_trace)
551 ! matrix_dcom has the negative sign and we include the negative sign of the coordinate
552 aat_tmp = aat_tmp + aat_prefactor*tmp_trace*lc_tmp*(-vcd_env%magnetic_origin_atom(gamma))
553 END DO
554 END DO
555
556 aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
557 = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
558 END DO
559
560 ! Nuclear contribution
561 CALL get_atomic_kind(particle_set(vcd_env%dcdr_env%lambda)%atomic_kind, kind_number=ikind)
562 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
563 IF (.NOT. ghost) THEN
564 DO alpha = 1, 3
565 aat_tmp = 0._dp
566 DO gamma = 1, 3
567 IF (levi_civita(alpha, gamma, vcd_env%dcdr_env%beta) == 0._dp) cycle
568 aat_tmp = aat_tmp + charge &
569 *levi_civita(alpha, gamma, vcd_env%dcdr_env%beta) &
570 *(particle_set(vcd_env%dcdr_env%lambda)%r(gamma) - vcd_env%magnetic_origin_atom(gamma))
571
572 aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
573 = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
574 END DO
575 END DO
576 END IF
577 END associate
578
579 CALL cp_fm_release(tmp_aomo)
580 CALL timestop(handle)
581 END SUBROUTINE aat_dv
582
583! **************************************************************************************************
584!> \brief Compute E_{alpha beta}^lambda = d/dV^lambda_beta <\mu_alpha> = d/dV^lambda_beta < \dot{r} >
585!> The directions alpha, beta are stored in vcd_env%dcdr_env
586!> \param vcd_env ...
587!> \param qs_env ...
588!> \author Edward Ditler, Tomas Zimmermann
589! **************************************************************************************************
590 SUBROUTINE apt_dv(vcd_env, qs_env)
591 TYPE(vcd_env_type) :: vcd_env
592 TYPE(qs_environment_type), POINTER :: qs_env
593
594 CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dV'
595 INTEGER, PARAMETER :: ispin = 1
596 REAL(dp), PARAMETER :: f_spin = 2._dp
597
598 INTEGER :: alpha, handle, ikind, nao, nmo
599 LOGICAL :: ghost
600 REAL(dp) :: charge
601 REAL(KIND=dp) :: apt_dcom, apt_difdip, apt_dipvel, &
602 apt_hcom, apt_rcom
603 TYPE(cp_fm_type) :: buf, matrix_dSdV_mo
604 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
605 POINTER :: sab_all
606 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
607 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
608
609 CALL timeset(routinen, handle)
610
611 CALL get_qs_env(qs_env=qs_env, &
612 sab_all=sab_all, &
613 particle_set=particle_set, &
614 qs_kind_set=qs_kind_set)
615
616 nmo = vcd_env%dcdr_env%nmo(ispin)
617 nao = vcd_env%dcdr_env%nao
618
619 associate(apt_el => vcd_env%apt_el_nvpt, &
620 apt_nuc => vcd_env%apt_nuc_nvpt, &
621 apt_total => vcd_env%apt_total_nvpt, &
622 mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), &
623 deltar => vcd_env%dcdr_env%deltaR)
624
625 ! build the full matrices
626 CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.true.)
627 CALL cp_fm_create(matrix_dsdv_mo, vcd_env%dcdr_env%momo_fm_struct(ispin)%struct)
628
629 ! STEP 1: dCV contribution (dipvel + commutator)
630 ! <mu|∂_alpha|nu> and <mu|[r_alpha, V]|nu> in AO basis
631 ! We compute tr(c_1^* x ∂_munu x c_0) + tr(c_0 x ∂_munu x c_1)
632 ! We compute tr(c_1^* x [,]_munu x c_0) + tr(c_0 x [,]_munu x c_1)
633 CALL cp_fm_scale_and_add(0._dp, vcd_env%dCV_prime(ispin), -1._dp, vcd_env%dCV(ispin))
634
635 ! Ref independent
636 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
637 buf, ncol=nmo)
638 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
639 1.0_dp, mo_coeff, buf, &
640 0.0_dp, matrix_dsdv_mo)
641
642 CALL parallel_gemm("N", "N", nao, nmo, nmo, &
643 -0.5_dp, mo_coeff, matrix_dsdv_mo, &
644 1.0_dp, vcd_env%dCV_prime(ispin))
645
646 ! + i∂ - i[Vnl, r]
647 DO alpha = 1, 3
648 CALL cp_fm_set_all(buf, 0.0_dp)
649 apt_dipvel = 0.0_dp
650
651 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao(alpha)%matrix, mo_coeff, buf, ncol=nmo)
652 CALL cp_fm_trace(buf, vcd_env%dCV_prime(ispin), apt_dipvel)
653 ! dipvel_ao = + < a | ∂ | b >
654 ! mo_coeff * dCV_prime = + iP1
655 apt_dipvel = 2._dp*apt_dipvel
656 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
657 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_dipvel
658 END DO
659
660 DO alpha = 1, 3
661 CALL cp_fm_set_all(buf, 0.0_dp)
662 apt_hcom = 0.0_dp
663 CALL cp_dbcsr_sm_fm_multiply(vcd_env%hcom(alpha)%matrix, mo_coeff, buf, ncol=nmo)
664 CALL cp_fm_trace(buf, vcd_env%dCV_prime(ispin), apt_hcom)
665
666 ! hcom = < a | [r, V] | b > = - < a | [V, r] | b >
667 ! mo_coeff * dCV_prime = + iP1
668 apt_hcom = +2._dp*apt_hcom
669
670 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
671 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_hcom
672 END DO !x/y/z
673
674 ! STEP 2: basis function derivative contribution
675 !! difdip_s
676 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, 0.0_dp)
677 CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s1(1)%matrix, &
678 vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix)
679 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, &
680 vcd_env%dcdr_env%lambda, direction_or=.true.)
681
682 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
683 buf, ncol=nmo, alpha=1._dp, beta=0._dp)
684 CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
685
686 apt_difdip = -f_spin*apt_difdip
687 apt_el(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) &
688 = apt_el(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) + apt_difdip
689
690 !! difdip(j, idir) = < a | r_j | ∂_idir b >
691 !! matrix_difdip2(beta, alpha) = < a | r_beta | ∂_alpha b >
692 ! matrix_difdip2 stores nuclear derivatives.
693 DO alpha = 1, 3 ! x/y/z for differentiated AO
694 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_difdip2(vcd_env%dcdr_env%beta, alpha)%matrix, mo_coeff, &
695 buf, ncol=nmo, alpha=1._dp, beta=0._dp)
696
697 CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
698 apt_difdip = -f_spin*apt_difdip
699 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
700 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + apt_difdip
701
702 END DO !alpha
703
704 ! STEP 3: The terms r * [V, r]
705 ! vcd_env%matrix_rrcom(alpha, beta) = r_beta * [V, r_alpha]
706 ! vcd_env%matrix_rcomr(alpha, beta) = [V, r_alpha] * r_beta
707 DO alpha = 1, 3 ! x/y/z for differentiated AO
708 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rcomr(alpha, vcd_env%dcdr_env%beta)%matrix)
709 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rrcom(alpha, vcd_env%dcdr_env%beta)%matrix)
710
711 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
712 sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
713 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
714 sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
715
716 CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
717 1.0_dp, -1.0_dp)
718
719 CALL cp_fm_set_all(buf, 0.0_dp)
720 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
721 CALL cp_fm_trace(mo_coeff, buf, apt_rcom)
722
723 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
724 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_rcom
725 END DO !alpha
726
727 ! STEP 4: pseudopotential derivative contribution
728 ! vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta) = - < mu | [ [V, r_beta], r_alpha ] | nu >
729 DO alpha = 1, 3 !x/y/z for differentiated AO
730 CALL cp_fm_set_all(buf, 0.0_dp)
731 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dcom(alpha, vcd_env%dcdr_env%beta)%matrix, mo_coeff, buf, ncol=nmo)
732 CALL cp_fm_trace(mo_coeff, buf, apt_dcom)
733 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
734 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_dcom
735 END DO !alpha
736
737 ! The reference point dependent terms:
738 !! difdip_munu
739 ! The additional term here is < a | db/dr(alpha)> * (delta_a - delta_b) * ref_point(beta)
740 ! in qs_env%matrix_s1(2:4) there is < da/dR | b > = - < da/dr | b > = < a | db/dr >
741 DO alpha = 1, 3
742 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, 0._dp)
743 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, 0._dp)
744 CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(alpha + 1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
745 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
746
747 ! < a | db/dr(alpha) > * R^lambda_beta * delta^lambda_nu
748 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
749 vcd_env%dcdr_env%lambda, direction_or=.true.)
750 ! < a | db/dr(alpha) > * R^lambda_beta * delta^lambda_mu
751 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
752 vcd_env%dcdr_env%lambda, direction_or=.false.)
753
754 ! < a | db/dr > * R^lambda_beta * ( delta^lambda_mu - delta^lambda_nu )
755 CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, &
756 1._dp, -1._dp)
757
758 CALL cp_fm_set_all(buf, 0.0_dp)
759 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, mo_coeff, buf, ncol=nmo)
760 CALL cp_fm_trace(mo_coeff, buf, apt_difdip)
761
762 ! And the whole contribution is
763 ! - < a | db/dr > * (mu - nu) * ref_point
764 apt_difdip = -apt_difdip*vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)
765
766 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
767 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_difdip
768 END DO
769
770 ! And the additional factor to rcom
771 ! < mu | [V, r] | nu > * R^lambda_beta * delta^lambda_mu
772 ! - < mu | [V, r] | nu > * R^lambda_beta * delta^lambda_nu
773 !
774 ! vcd_env%hcom(alpha) = - < mu | [V, r_alpha] | nu >
775 ! particle_set(lambda)%r(vcd_env%dcdr_env%beta) = R^lambda_beta
776
777 DO alpha = 1, 3
778 CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, 0._dp)
779 CALL dbcsr_desymmetrize(vcd_env%hcom(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
780 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix)
781
782 ! < mu | [V, r] | nu > * delta^lambda_nu
783 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
784 vcd_env%dcdr_env%lambda, direction_or=.true.)
785 ! < mu | [V, r] | nu > * delta^lambda_mu
786 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, qs_kind_set, "ORB", sab_all, &
787 vcd_env%dcdr_env%lambda, direction_or=.false.)
788
789 CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, &
790 vcd_env%dcdr_env%matrix_nosym_temp2(alpha)%matrix, -1._dp, +1._dp)
791
792 CALL cp_fm_set_all(buf, 0.0_dp)
793 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(alpha)%matrix, mo_coeff, buf, ncol=nmo)
794 CALL cp_fm_trace(mo_coeff, buf, apt_rcom)
795 apt_rcom = -vcd_env%spatial_origin_atom(vcd_env%dcdr_env%beta)*apt_rcom
796
797 apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
798 = apt_el(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + f_spin*apt_rcom
799 END DO
800
801 ! STEP 5: nuclear contribution
802 associate(atomic_kind => particle_set(vcd_env%dcdr_env%lambda)%atomic_kind)
803 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
804 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
805 IF (.NOT. ghost) THEN
806 apt_nuc(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) = &
807 apt_nuc(vcd_env%dcdr_env%beta, vcd_env%dcdr_env%beta, vcd_env%dcdr_env%lambda) + charge
808 END IF
809 END associate
810
811 ! STEP 6: deallocations
812 CALL cp_fm_release(buf)
813 CALL cp_fm_release(matrix_dsdv_mo)
814
815 END associate
816
817 CALL timestop(handle)
818 END SUBROUTINE apt_dv
819
820! **************************************************************************************************
821!> \brief Initialize the matrices for the NVPT calculation
822!> \param vcd_env ...
823!> \param qs_env ...
824!> \author Edward Ditler
825! **************************************************************************************************
826 SUBROUTINE prepare_per_atom_vcd(vcd_env, qs_env)
827 TYPE(vcd_env_type) :: vcd_env
828 TYPE(qs_environment_type), POINTER :: qs_env
829
830 CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_per_atom_vcd'
831
832 INTEGER :: handle, i, ispin, j
833 TYPE(cell_type), POINTER :: cell
834 TYPE(dft_control_type), POINTER :: dft_control
835 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
836 POINTER :: sab_all, sab_orb, sap_ppnl
837 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
838 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
839
840 CALL timeset(routinen, handle)
841
842 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
843 sab_orb=sab_orb, sab_all=sab_all, sap_ppnl=sap_ppnl, &
844 qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
845
846 IF (vcd_env%distributed_origin) THEN
847 vcd_env%magnetic_origin_atom(:) = particle_set(vcd_env%dcdr_env%lambda)%r(:) - vcd_env%magnetic_origin(:)
848 vcd_env%spatial_origin_atom = particle_set(vcd_env%dcdr_env%lambda)%r(:) - vcd_env%spatial_origin(:)
849 END IF
850
851 ! Reset the matrices
852 DO ispin = 1, dft_control%nspins
853 DO j = 1, 3
854 CALL dbcsr_set(vcd_env%matrix_dSdV(j)%matrix, 0._dp)
855 CALL dbcsr_set(vcd_env%matrix_drpnl(j)%matrix, 0._dp)
856
857 DO i = 1, 3
858 CALL dbcsr_set(vcd_env%matrix_dcom(i, j)%matrix, 0.0_dp)
859 CALL dbcsr_set(vcd_env%matrix_difdip2(i, j)%matrix, 0._dp)
860 END DO
861 END DO
862 CALL cp_fm_set_all(vcd_env%op_dV(ispin), 0._dp)
863 CALL dbcsr_set(vcd_env%matrix_hxc_dsdv(ispin)%matrix, 0._dp)
864 END DO
865
866 ! operator dV
867 ! <mu|d/dV_beta [V, r_alpha]|nu>
868 CALL build_dcom_rpnl(vcd_env%matrix_dcom, qs_kind_set, sab_orb, sap_ppnl, &
869 dft_control%qs_control%eps_ppnl, particle_set, vcd_env%dcdr_env%lambda)
870
871 ! PP derivative. build_com_mom_nl returns [r, Vnl], while matrix_drpnl
872 ! historically stores [Vnl, r] = -[r, Vnl].
873 CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
874 particle_set, cell=cell, matrix_rv=vcd_env%matrix_drpnl, &
875 pseudoatom=vcd_env%dcdr_env%lambda)
876 DO j = 1, 3
877 CALL dbcsr_scale(vcd_env%matrix_drpnl(j)%matrix, alpha_scalar=-1._dp)
878 END DO
879 ! lin_mom
880 DO i = 1, 3
881 CALL dbcsr_set(vcd_env%dipvel_ao_delta(i)%matrix, 0._dp)
882 CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dipvel_ao(i)%matrix)
883 END DO
884
885 CALL hr_mult_by_delta_3d(vcd_env%dipvel_ao_delta, qs_kind_set, "ORB", &
886 sab_all, vcd_env%dcdr_env%delta_basis_function, direction_or=.true.)
887
888 ! dS/dV
889 CALL build_dsdv_matrix(qs_env, vcd_env%matrix_dSdV, &
890 deltar=vcd_env%dcdr_env%delta_basis_function, &
891 rcc=vcd_env%spatial_origin_atom)
892
893 CALL build_local_moments_der_matrix(qs_env, vcd_env%matrix_difdip2, 1, 0, &
894 ref_point=[0._dp, 0._dp, 0._dp], basis_type="ORB", &
895 ordered=.true., lambda=vcd_env%dcdr_env%lambda)
896 ! AAT
897 ! moments_throw: x, y, z, xx, xy, xz, yy, yz, zz
898 ! moments_der: (moment, xyz derivative)
899 ! build_local_moments_der_matrix uses adbdr for calculating derivatives of the *primitive*
900 ! on the right. So the resulting
901 ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
902 DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
903 DO j = 1, 3
904 CALL dbcsr_set(vcd_env%moments_der_right(i, j)%matrix, 0.0_dp)
905 CALL dbcsr_set(vcd_env%moments_der_left(i, j)%matrix, 0.0_dp)
906 END DO
907 END DO
908
909 DO i = 1, 9
910 DO j = 1, 3 ! derivatives
911 CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_right(i, j)%matrix) ! A2
912 CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_left(i, j)%matrix) ! A1
913
914 ! - < mu | r_beta r_gamma ∂_delta | nu > * (mu/nu == lambda)
915 CALL hr_mult_by_delta_1d(vcd_env%moments_der_right(i, j)%matrix, qs_kind_set, "ORB", &
916 sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
917 CALL hr_mult_by_delta_1d(vcd_env%moments_der_left(i, j)%matrix, qs_kind_set, "ORB", &
918 sab_all, direction_or=.false., lambda=vcd_env%dcdr_env%lambda)
919 END DO
920 END DO
921
922 DO i = 1, 3
923 DO j = 1, 3
924 CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
925 END DO
926 END DO
927
928 CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
929 particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
930 matrix_r_doublecom=vcd_env%matrix_r_doublecom, &
931 pseudoatom=vcd_env%dcdr_env%lambda)
932
933 CALL timestop(handle)
934
935 END SUBROUTINE prepare_per_atom_vcd
936
937! **************************************************************************************************
938!> \brief What we are building here is the operator for the NVPT response:
939!> H0 * C1 - S0 * E0 * C1 = - op_dV
940!> linres_solver = - [ H1 * C0 - S1 * C0 * E0 ]
941!> with
942!> H1 * C0 = dH/dV * C0
943!> + i[∂]δ * C0
944!> - i S0 * C^(1,R)
945!> + i S0 * C0 * (C0 * S^(1,R) * C0)
946!> - S1 * C0 * E0
947!>
948!> H1 * C0 = + i (Hr - rH) * C0 [STEP 1]
949!> + i[∂]δ * C0 [STEP 2]
950!> - i[V, r]δ * C0 [STEP 3]
951!> - i S0 * C^(1,R) [STEP 4]
952!> - S1 * C0 * E0 [STEP 5]
953!> \param vcd_env ...
954!> \param qs_env ...
955!> \author Edward Ditler, Tomas Zimmermann
956! **************************************************************************************************
957 SUBROUTINE vcd_build_op_dv(vcd_env, qs_env)
958 TYPE(vcd_env_type) :: vcd_env
959 TYPE(qs_environment_type), POINTER :: qs_env
960
961 CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_build_op_dV'
962 INTEGER, PARAMETER :: ispin = 1
963
964 INTEGER :: handle, nao, nmo
965 TYPE(cp_fm_type) :: buf
966 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
967 POINTER :: sab_all
968 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
969
970 CALL timeset(routinen, handle)
971
972 CALL get_qs_env(qs_env=qs_env, &
973 sab_all=sab_all, &
974 qs_kind_set=qs_kind_set)
975
976 nmo = vcd_env%dcdr_env%nmo(1)
977 nao = vcd_env%dcdr_env%nao
978
979 CALL build_matrix_hr_rh(vcd_env, qs_env, vcd_env%spatial_origin_atom)
980
981 ! STEP 1: hr-rh
982 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_hr(ispin, vcd_env%dcdr_env%beta)%matrix)
983 CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rh(ispin, vcd_env%dcdr_env%beta)%matrix)
984
985 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
986 sab_all, vcd_env%dcdr_env%lambda, direction_or=.true.)
987 CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
988 sab_all, vcd_env%dcdr_env%lambda, direction_or=.false.)
989 CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
990 vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
991 1.0_dp, -1.0_dp)
992
993 associate(mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
994 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, &
995 vcd_env%op_dV(ispin), ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
996
997 ! STEP 2: electronic momentum operator contribution
998 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao_delta(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
999 vcd_env%op_dV(ispin), &
1000 ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
1001
1002 ! STEP 3: +dV_ppnl/dV, but matrix_drpnl stores the negative of dV_ppnl
1003 ! The arguments (-1, 1) are swapped wrt to the hr-rh term, implying that
1004 ! direction_Or and direction_hr do what they should.
1005 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_drpnl(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
1006 vcd_env%op_dV(ispin), &
1007 ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
1008
1009 ! STEP 4: - S0 * C^(1,R)
1010 CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1011 CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), &
1012 vcd_env%op_dV(1), ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
1013
1014 ! STEP 5: -S(1,V) * C0 * E0
1015 CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
1016 buf, nmo, alpha=1.0_dp, beta=0.0_dp)
1017 CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
1018 -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &
1019 1.0_dp, vcd_env%op_dV(ispin))
1020
1021 CALL cp_fm_release(buf)
1022 END associate
1023
1024 ! We have built op_dV but plug -op_dV into the linres_solver
1025 CALL cp_fm_scale(-1.0_dp, vcd_env%op_dV(1))
1026
1027 ! Revert the matrices
1028 CALL build_matrix_hr_rh(vcd_env, qs_env, [0._dp, 0._dp, 0._dp])
1029
1030 CALL timestop(handle)
1031 END SUBROUTINE vcd_build_op_dv
1032
1033! *****************************************************************************
1034!> \brief Get the dC/dV using the vcd_env%op_dV
1035!> \param vcd_env ...
1036!> \param p_env ...
1037!> \param qs_env ...
1038!> \author Edward Ditler, Tomas Zimmermann
1039! **************************************************************************************************
1040 SUBROUTINE vcd_response_dv(vcd_env, p_env, qs_env)
1041
1042 TYPE(vcd_env_type) :: vcd_env
1043 TYPE(qs_p_env_type) :: p_env
1044 TYPE(qs_environment_type), POINTER :: qs_env
1045
1046 CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_response_dV'
1047 INTEGER, PARAMETER :: ispin = 1
1048
1049 INTEGER :: handle, output_unit
1050 LOGICAL :: failure, should_stop
1051 TYPE(cp_fm_type), DIMENSION(1) :: h1_psi0, psi1
1052 TYPE(cp_logger_type), POINTER :: logger
1053 TYPE(linres_control_type), POINTER :: linres_control
1054 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1055 TYPE(section_vals_type), POINTER :: lr_section, vcd_section
1056
1057 CALL timeset(routinen, handle)
1058 failure = .false.
1059
1060 NULLIFY (linres_control, lr_section, logger)
1061
1062 CALL get_qs_env(qs_env=qs_env, &
1063 linres_control=linres_control, &
1064 mos=mos)
1065
1066 logger => cp_get_default_logger()
1067 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
1068 vcd_section => section_vals_get_subs_vals(qs_env%input, &
1069 "PROPERTIES%LINRES%vcd")
1070
1071 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
1072 extension=".linresLog")
1073 IF (output_unit > 0) THEN
1074 WRITE (unit=output_unit, fmt="(T10,A,/)") &
1075 "*** Self consistent optimization of the response wavefunction ***"
1076 END IF
1077
1078 associate(psi0_order => vcd_env%dcdr_env%mo_coeff)
1079 CALL cp_fm_create(psi1(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.true.)
1080 CALL cp_fm_create(h1_psi0(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1081
1082 ! Restart
1083 IF (linres_control%linres_restart) THEN
1084 CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
1085 ELSE
1086 CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
1087 END IF
1088
1089 IF (output_unit > 0) THEN
1090 WRITE (output_unit, *) &
1091 "Response to the perturbation operator referring to the velocity of atom ", &
1092 vcd_env%dcdr_env%lambda, " in "//achar(vcd_env%dcdr_env%beta + 119)
1093 END IF
1094
1095 ! First response to get dCR
1096 ! (H0-E0) psi1 = (H1-E1) psi0
1097 ! psi1 = the perturbed wavefunction
1098 ! h1_psi0 = (H1-E1)
1099 ! psi0_order = the unperturbed wavefunction
1100 ! Second response to get dCV
1101 CALL cp_fm_set_all(vcd_env%dCV(ispin), 0.0_dp)
1102 CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
1103 CALL cp_fm_to_fm(vcd_env%op_dV(ispin), h1_psi0(ispin))
1104
1105 linres_control%lr_triplet = .false. ! we do singlet response
1106 linres_control%do_kernel = .false. ! no coupled response since imaginary perturbation
1107 linres_control%converged = .false.
1108 CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
1109 output_unit, should_stop)
1110 CALL cp_fm_to_fm(psi1(ispin), vcd_env%dCV(ispin))
1111
1112 ! Write the new result to the restart file
1113 IF (linres_control%linres_restart) THEN
1114 CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
1115 END IF
1116
1117 END associate
1118
1119 ! clean up
1120 CALL cp_fm_release(psi1(ispin))
1121 CALL cp_fm_release(h1_psi0(ispin))
1122
1123 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1124 "PRINT%PROGRAM_RUN_INFO")
1125
1126 CALL timestop(handle)
1127 END SUBROUTINE vcd_response_dv
1128
1129END MODULE qs_vcd
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:78
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
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:76
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_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
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:799
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:170
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:1502
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:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
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:1073
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:735
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:631
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, mimic, 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, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
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, cneo_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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, 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:491
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, silent)
scf loop to optimize the first order wavefunctions (psi1) given a perturbation as an operator applied...
Definition qs_linres_methods.F:218
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::build_local_moments_der_matrix
subroutine, public build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, ref_point, moments, basis_type, minimum_image, ordered, lambda, deltar, neighbor_image)
Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b > Optionally sto...
Definition qs_moments.F:487
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:1725
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:2260
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:1426
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:117
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:635
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:485
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:780
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
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:224
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
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:53
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