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