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