Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : MODULE qs_vcd
8 : USE atomic_kind_types, ONLY: get_atomic_kind
9 : USE cell_types, ONLY: cell_type
10 : USE commutator_rpnl, ONLY: build_com_mom_nl
11 : USE cp_control_types, ONLY: dft_control_type
12 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
13 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
14 : cp_fm_scale_and_add,&
15 : cp_fm_trace
16 : USE cp_fm_types, ONLY: cp_fm_create,&
17 : cp_fm_release,&
18 : cp_fm_set_all,&
19 : cp_fm_to_fm,&
20 : cp_fm_type
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type
23 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
24 : cp_print_key_unit_nr
25 : USE dbcsr_api, ONLY: dbcsr_add,&
26 : dbcsr_copy,&
27 : dbcsr_desymmetrize,&
28 : dbcsr_set
29 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
30 : section_vals_type
31 : USE kinds, ONLY: dp
32 : USE parallel_gemm_api, ONLY: parallel_gemm
33 : USE particle_types, ONLY: particle_type
34 : USE qs_dcdr_ao, ONLY: hr_mult_by_delta_1d
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_kind_types, ONLY: get_qs_kind,&
38 : qs_kind_type
39 : USE qs_linres_methods, ONLY: linres_solver
40 : USE qs_linres_types, ONLY: linres_control_type,&
41 : vcd_env_type
42 : USE qs_mo_types, ONLY: mo_set_type
43 : USE qs_moments, ONLY: dipole_velocity_deriv
44 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
45 : USE qs_p_env_types, ONLY: qs_p_env_type
46 : USE qs_vcd_ao, ONLY: build_dSdV_matrix,&
47 : build_dcom_rpnl,&
48 : build_drpnl_matrix,&
49 : build_matrix_hr_rh,&
50 : hr_mult_by_delta_3d
51 : USE qs_vcd_utils, ONLY: vcd_read_restart,&
52 : vcd_write_restart
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 : PUBLIC :: prepare_per_atom_vcd
59 : PUBLIC :: vcd_build_op_dV
60 : PUBLIC :: vcd_response_dV
61 : PUBLIC :: apt_dV
62 : PUBLIC :: aat_dV
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vcd'
65 :
66 : REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = RESHAPE((/ &
67 : 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
68 : 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
69 : 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
70 : (/3, 3, 3/))
71 : INTEGER, DIMENSION(3, 3), PARAMETER :: multipole_2d_to_1d = RESHAPE([4, 5, 6, 5, 7, 8, 6, 8, 9], [3, 3])
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Compute I_{alpha beta}^lambda = d/dV^lambda_beta <m_alpha> = d/dV^lambda_beta < r x \dot{r} >
76 : !> The directions alpha, beta are stored in vcd_env%dcdr_env
77 : !> \param vcd_env ...
78 : !> \param qs_env ...
79 : !> \author Edward Ditler
80 : ! **************************************************************************************************
81 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 .EQ. 1) aat_prefactor = aat_prefactor*2.0_dp
121 :
122 : ! The non-PP part of the AAT consists of four contributions:
123 : ! (A1): + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda)
124 : ! (A2): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta ∂_delta | nu > * (nu == lambda)
125 : ! (B): - P^0 * ε_{alpha gamma delta} * < mu | r_gamma | nu > * (delta == beta) * (nu == lambda)
126 : ! (C): + iP^1 * ε_{alpha gamma delta} * < mu | r_gamma ∂_delta | nu >
127 :
128 : ! (A1) + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda)
129 : ! (A2) - P^0 * ε_{alpha gamma delta} * < mu | r_gamma r_beta ∂_delta | nu > * (nu == lambda)
130 : ! Conjecture : It doesn't matter that the beta and gamma are swapped around!
131 : ! We define o = | ∂_delta nu >
132 : ! and then < a | r_beta r_gamma | o > = < a | r_gamma r_beta | o>
133 : ! (A) + P^0 * ε_{alpha gamma delta} * < mu | r_beta r_gamma ∂_delta | nu > * (mu == lambda - nu == lambda)
134 : ! We have built the matrices - < mu | r_beta r_gamma ∂_delta | nu > in vcd_env%moments_der
135 : ! moments_der(1:9; 1:3) = moments_der(x, y, z, xx, xy, xz, yy, yz, zz;
136 : ! x, y, z)
137 :
138 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)
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
875 : CALL build_drpnl_matrix(vcd_env%matrix_drpnl, qs_kind_set, sab_all, sap_ppnl, &
876 6 : dft_control%qs_control%eps_ppnl, particle_set, pseudoatom=vcd_env%dcdr_env%lambda)
877 : ! lin_mom
878 24 : DO i = 1, 3
879 18 : CALL dbcsr_set(vcd_env%dipvel_ao_delta(i)%matrix, 0._dp)
880 24 : CALL dbcsr_copy(vcd_env%dipvel_ao_delta(i)%matrix, vcd_env%dipvel_ao(i)%matrix)
881 : END DO
882 :
883 : CALL hr_mult_by_delta_3d(vcd_env%dipvel_ao_delta, qs_kind_set, "ORB", &
884 6 : sab_all, vcd_env%dcdr_env%delta_basis_function, direction_Or=.TRUE.)
885 :
886 : ! dS/dV
887 : CALL build_dSdV_matrix(qs_env, vcd_env%matrix_dSdV, &
888 : deltaR=vcd_env%dcdr_env%delta_basis_function, &
889 6 : rcc=vcd_env%spatial_origin_atom)
890 :
891 : CALL dipole_velocity_deriv(qs_env, vcd_env%matrix_difdip2, 1, lambda=vcd_env%dcdr_env%lambda, &
892 6 : rc=[0._dp, 0._dp, 0._dp])
893 : ! AAT
894 : ! moments_throw: x, y, z, xx, xy, xz, yy, yz, zz
895 : ! moments_der: (moment, xyz derivative)
896 : ! build_local_moments_der_matrix uses adbdr for calculating derivatives of the *primitive*
897 : ! on the right. So the resulting
898 : ! moments_der(moment, delta) = - < a | moment \partial_\delta | b >
899 60 : DO i = 1, 9 ! x, y, z, xx, xy, xz, yy, yz, zz
900 222 : DO j = 1, 3
901 162 : CALL dbcsr_set(vcd_env%moments_der_right(i, j)%matrix, 0.0_dp)
902 216 : CALL dbcsr_set(vcd_env%moments_der_left(i, j)%matrix, 0.0_dp)
903 : END DO
904 : END DO
905 :
906 60 : DO i = 1, 9
907 222 : DO j = 1, 3 ! derivatives
908 162 : CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_right(i, j)%matrix) ! A2
909 162 : CALL dbcsr_desymmetrize(vcd_env%moments_der(i, j)%matrix, vcd_env%moments_der_left(i, j)%matrix) ! A1
910 :
911 : ! - < mu | r_beta r_gamma ∂_delta | nu > * (mu/nu == lambda)
912 : CALL hr_mult_by_delta_1d(vcd_env%moments_der_right(i, j)%matrix, qs_kind_set, "ORB", &
913 162 : sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
914 : CALL hr_mult_by_delta_1d(vcd_env%moments_der_left(i, j)%matrix, qs_kind_set, "ORB", &
915 216 : sab_all, direction_Or=.FALSE., lambda=vcd_env%dcdr_env%lambda)
916 : END DO
917 : END DO
918 :
919 24 : DO i = 1, 3
920 78 : DO j = 1, 3
921 72 : CALL dbcsr_set(vcd_env%matrix_r_doublecom(i, j)%matrix, 0._dp)
922 : END DO
923 : END DO
924 :
925 : CALL build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, dft_control%qs_control%eps_ppnl, &
926 : particle_set, ref_point=[0._dp, 0._dp, 0._dp], cell=cell, &
927 : matrix_r_doublecom=vcd_env%matrix_r_doublecom, &
928 6 : pseudoatom=vcd_env%dcdr_env%lambda)
929 :
930 6 : CALL timestop(handle)
931 :
932 6 : END SUBROUTINE prepare_per_atom_vcd
933 :
934 : ! **************************************************************************************************
935 : !> \brief What we are building here is the operator for the NVPT response:
936 : !> H0 * C1 - S0 * E0 * C1 = - op_dV
937 : !> linres_solver = - [ H1 * C0 - S1 * C0 * E0 ]
938 : !> with
939 : !> H1 * C0 = dH/dV * C0
940 : !> + i[∂]δ * C0
941 : !> - i S0 * C^(1,R)
942 : !> + i S0 * C0 * (C0 * S^(1,R) * C0)
943 : !> - S1 * C0 * E0
944 : !>
945 : !> H1 * C0 = + i (Hr - rH) * C0 [STEP 1]
946 : !> + i[∂]δ * C0 [STEP 2]
947 : !> - i[V, r]δ * C0 [STEP 3]
948 : !> - i S0 * C^(1,R) [STEP 4]
949 : !> - S1 * C0 * E0 [STEP 5]
950 : !> \param vcd_env ...
951 : !> \param qs_env ...
952 : !> \author Edward Ditler, Tomas Zimmermann
953 : ! **************************************************************************************************
954 18 : SUBROUTINE vcd_build_op_dV(vcd_env, qs_env)
955 : TYPE(vcd_env_type) :: vcd_env
956 : TYPE(qs_environment_type), POINTER :: qs_env
957 :
958 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_build_op_dV'
959 : INTEGER, PARAMETER :: ispin = 1
960 :
961 : INTEGER :: handle, nao, nmo
962 : TYPE(cp_fm_type) :: buf
963 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
964 18 : POINTER :: sab_all
965 18 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
966 :
967 18 : CALL timeset(routineN, handle)
968 :
969 : CALL get_qs_env(qs_env=qs_env, &
970 : sab_all=sab_all, &
971 18 : qs_kind_set=qs_kind_set)
972 :
973 18 : nmo = vcd_env%dcdr_env%nmo(1)
974 18 : nao = vcd_env%dcdr_env%nao
975 :
976 18 : CALL build_matrix_hr_rh(vcd_env, qs_env, vcd_env%spatial_origin_atom)
977 :
978 : ! STEP 1: hr-rh
979 18 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_hr(ispin, vcd_env%dcdr_env%beta)%matrix)
980 18 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rh(ispin, vcd_env%dcdr_env%beta)%matrix)
981 :
982 : CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
983 18 : sab_all, vcd_env%dcdr_env%lambda, direction_or=.TRUE.)
984 : CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, qs_kind_set, "ORB", &
985 18 : sab_all, vcd_env%dcdr_env%lambda, direction_or=.FALSE.)
986 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
987 : vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
988 18 : 1.0_dp, -1.0_dp)
989 :
990 : ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
991 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, &
992 18 : vcd_env%op_dV(ispin), ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
993 :
994 : ! STEP 2: electronic momentum operator contribution
995 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dipvel_ao_delta(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
996 : vcd_env%op_dV(ispin), &
997 18 : ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
998 :
999 : ! STEP 3: +dV_ppnl/dV, but build_drpnl_matrix gives the negative of dV_ppnl
1000 : ! The arguments (-1, 1) are swapped wrt to the hr-rh term, implying that
1001 : ! direction_Or and direction_hr do what they should.
1002 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_drpnl(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
1003 : vcd_env%op_dV(ispin), &
1004 18 : ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
1005 :
1006 : ! STEP 4: - S0 * C^(1,R)
1007 18 : CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1008 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s1(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), &
1009 18 : vcd_env%op_dV(1), ncol=nmo, alpha=-1.0_dp, beta=1.0_dp)
1010 :
1011 : ! STEP 5: -S(1,V) * C0 * E0
1012 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdV(vcd_env%dcdr_env%beta)%matrix, mo_coeff, &
1013 18 : buf, nmo, alpha=1.0_dp, beta=0.0_dp)
1014 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
1015 : -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), &
1016 18 : 1.0_dp, vcd_env%op_dV(ispin))
1017 :
1018 36 : CALL cp_fm_release(buf)
1019 : END ASSOCIATE
1020 :
1021 : ! We have built op_dV but plug -op_dV into the linres_solver
1022 18 : CALL cp_fm_scale(-1.0_dp, vcd_env%op_dV(1))
1023 :
1024 : ! Revert the matrices
1025 18 : CALL build_matrix_hr_rh(vcd_env, qs_env, [0._dp, 0._dp, 0._dp])
1026 :
1027 18 : CALL timestop(handle)
1028 36 : END SUBROUTINE vcd_build_op_dV
1029 :
1030 : ! *****************************************************************************
1031 : !> \brief Get the dC/dV using the vcd_env%op_dV
1032 : !> \param vcd_env ...
1033 : !> \param p_env ...
1034 : !> \param qs_env ...
1035 : !> \author Edward Ditler, Tomas Zimmermann
1036 : ! **************************************************************************************************
1037 18 : SUBROUTINE vcd_response_dV(vcd_env, p_env, qs_env)
1038 :
1039 : TYPE(vcd_env_type) :: vcd_env
1040 : TYPE(qs_p_env_type) :: p_env
1041 : TYPE(qs_environment_type), POINTER :: qs_env
1042 :
1043 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vcd_response_dV'
1044 : INTEGER, PARAMETER :: ispin = 1
1045 :
1046 : INTEGER :: handle, output_unit
1047 : LOGICAL :: failure, should_stop
1048 72 : TYPE(cp_fm_type), DIMENSION(1) :: h1_psi0, psi1
1049 : TYPE(cp_logger_type), POINTER :: logger
1050 : TYPE(linres_control_type), POINTER :: linres_control
1051 18 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1052 : TYPE(section_vals_type), POINTER :: lr_section, vcd_section
1053 :
1054 18 : CALL timeset(routineN, handle)
1055 18 : failure = .FALSE.
1056 :
1057 18 : NULLIFY (linres_control, lr_section, logger)
1058 :
1059 : CALL get_qs_env(qs_env=qs_env, &
1060 : linres_control=linres_control, &
1061 18 : mos=mos)
1062 :
1063 18 : logger => cp_get_default_logger()
1064 18 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
1065 : vcd_section => section_vals_get_subs_vals(qs_env%input, &
1066 18 : "PROPERTIES%LINRES%vcd")
1067 :
1068 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
1069 18 : extension=".linresLog")
1070 18 : IF (output_unit > 0) THEN
1071 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
1072 9 : "*** Self consistent optimization of the response wavefunction ***"
1073 : END IF
1074 :
1075 : ASSOCIATE (psi0_order => vcd_env%dcdr_env%mo_coeff)
1076 18 : CALL cp_fm_create(psi1(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1077 18 : CALL cp_fm_create(h1_psi0(ispin), vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
1078 :
1079 : ! Restart
1080 18 : IF (linres_control%linres_restart) THEN
1081 18 : CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
1082 : ELSE
1083 0 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
1084 : END IF
1085 :
1086 18 : IF (output_unit > 0) THEN
1087 : WRITE (output_unit, *) &
1088 9 : "Response to the perturbation operator referring to the velocity of atom ", &
1089 18 : vcd_env%dcdr_env%lambda, " in "//ACHAR(vcd_env%dcdr_env%beta + 119)
1090 : END IF
1091 :
1092 : ! First response to get dCR
1093 : ! (H0-E0) psi1 = (H1-E1) psi0
1094 : ! psi1 = the perturbed wavefunction
1095 : ! h1_psi0 = (H1-E1)
1096 : ! psi0_order = the unperturbed wavefunction
1097 : ! Second response to get dCV
1098 18 : CALL cp_fm_set_all(vcd_env%dCV(ispin), 0.0_dp)
1099 18 : CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
1100 18 : CALL cp_fm_to_fm(vcd_env%op_dV(ispin), h1_psi0(ispin))
1101 :
1102 18 : linres_control%lr_triplet = .FALSE. ! we do singlet response
1103 18 : linres_control%do_kernel = .FALSE. ! no coupled response since imaginary perturbation
1104 18 : linres_control%converged = .FALSE.
1105 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
1106 18 : output_unit, should_stop)
1107 18 : CALL cp_fm_to_fm(psi1(ispin), vcd_env%dCV(ispin))
1108 :
1109 : ! Write the new result to the restart file
1110 36 : IF (linres_control%linres_restart) THEN
1111 18 : CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, vcd_env%dcdr_env%beta, "dCdV")
1112 : END IF
1113 :
1114 : END ASSOCIATE
1115 :
1116 : ! clean up
1117 18 : CALL cp_fm_release(psi1(ispin))
1118 18 : CALL cp_fm_release(h1_psi0(ispin))
1119 :
1120 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1121 18 : "PRINT%PROGRAM_RUN_INFO")
1122 :
1123 18 : CALL timestop(handle)
1124 36 : END SUBROUTINE vcd_response_dV
1125 :
1126 : END MODULE qs_vcd
|
