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_mfp
8 : USE atomic_kind_types, ONLY: get_atomic_kind
9 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
10 : gto_basis_set_type
11 : USE cell_types, ONLY: cell_type
12 : USE commutator_rpnl, ONLY: build_com_vnl_giao
13 : USE cp_control_types, ONLY: dft_control_type
14 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
15 : dbcsr_copy,&
16 : dbcsr_desymmetrize,&
17 : dbcsr_get_block_p,&
18 : dbcsr_p_type,&
19 : dbcsr_scale,&
20 : dbcsr_set,&
21 : dbcsr_type
22 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
23 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_release,&
26 : cp_fm_set_all,&
27 : cp_fm_to_fm,&
28 : cp_fm_type
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_type
31 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
32 : cp_print_key_unit_nr
33 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
34 : section_vals_type
35 : USE kinds, ONLY: dp
36 : USE parallel_gemm_api, ONLY: parallel_gemm
37 : USE particle_types, ONLY: particle_type
38 : USE qs_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : USE qs_integral_utils, ONLY: basis_set_list_setup
41 : USE qs_kind_types, ONLY: get_qs_kind,&
42 : qs_kind_type
43 : USE qs_linres_methods, ONLY: linres_solver
44 : USE qs_linres_types, ONLY: linres_control_type,&
45 : vcd_env_type
46 : USE qs_mo_types, ONLY: mo_set_type
47 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
48 : get_neighbor_list_set_p,&
49 : neighbor_list_iterate,&
50 : neighbor_list_iterator_create,&
51 : neighbor_list_iterator_p_type,&
52 : neighbor_list_iterator_release,&
53 : neighbor_list_set_p_type
54 : USE qs_p_env_types, ONLY: qs_p_env_type
55 : USE qs_vcd_ao, ONLY: hr_mult_by_delta_1d
56 : USE qs_vcd_utils, ONLY: vcd_read_restart,&
57 : vcd_write_restart
58 :
59 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
60 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
61 : !$ omp_init_lock, omp_set_lock, &
62 : !$ omp_unset_lock, omp_destroy_lock
63 : #include "./base/base_uses.f90"
64 :
65 : IMPLICIT NONE
66 :
67 : PRIVATE
68 : PUBLIC :: mfp_build_operator_gauge_independent
69 : PUBLIC :: mfp_build_operator_gauge_dependent
70 : PUBLIC :: mfp_response
71 : PUBLIC :: mfp_aat
72 : PUBLIC :: multiply_by_position
73 :
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mfp'
75 :
76 : REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = RESHAPE((/ &
77 : 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, &
78 : 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, &
79 : 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/), &
80 : (/3, 3, 3/))
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief ...
85 : !> \param vcd_env ...
86 : !> \param qs_env ...
87 : !> \author Edward Ditler
88 : ! **************************************************************************************************
89 0 : SUBROUTINE mfp_aat(vcd_env, qs_env)
90 : TYPE(vcd_env_type) :: vcd_env
91 : TYPE(qs_environment_type), POINTER :: qs_env
92 :
93 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mfp_aat'
94 : INTEGER, PARAMETER :: ispin = 1
95 : REAL(dp), DIMENSION(3), PARAMETER :: gauge_origin = 0._dp
96 : REAL(dp), PARAMETER :: f_spin = 2._dp
97 :
98 : INTEGER :: alpha, delta, gamma, handle, ikind, nao, &
99 : nmo
100 : LOGICAL :: ghost
101 : REAL(dp) :: aat_linmom, aat_moment, aat_moment_der, &
102 : aat_overlap, aat_tmp, charge, lc_tmp, &
103 : tmp_trace
104 : TYPE(cp_fm_type) :: buf, buf2, matrix_dSdB_mo
105 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
106 0 : POINTER :: sab_all
107 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
109 :
110 0 : CALL timeset(routineN, handle)
111 :
112 : ! Some setup
113 0 : nmo = vcd_env%dcdr_env%nmo(ispin)
114 0 : nao = vcd_env%dcdr_env%nao
115 :
116 : ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin), aat_atom => vcd_env%aat_atom_mfp)
117 :
118 : CALL get_qs_env(qs_env=qs_env, &
119 : sab_all=sab_all, &
120 : qs_kind_set=qs_kind_set, &
121 0 : particle_set=particle_set)
122 :
123 : ! Set up the matrices
124 0 : CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
125 0 : CALL cp_fm_create(buf2, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
126 0 : CALL cp_fm_create(matrix_dSdB_mo, vcd_env%dcdr_env%momo_fm_struct(ispin)%struct)
127 :
128 : ! First prepare the coefficients dCB_prime = -dCB - 0.5 S1_ij
129 : ! The first negative because we get the negative coefficients out of the linres solver
130 : ! The second term due to the occ-occ block of the U matrix in C1 = U * C0
131 0 : DO alpha = 1, 3
132 : ! CALL cp_fm_scale_and_add(0._dp, vcd_env%dCB_prime(alpha), -1._dp, vcd_env%dCB(alpha))
133 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdB(alpha)%matrix, mo_coeff, &
134 0 : buf, ncol=nmo, alpha=1._dp, beta=0._dp)
135 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
136 : 1.0_dp, mo_coeff, buf, &
137 0 : 0.0_dp, matrix_dSdB_mo)
138 :
139 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
140 : -0.5_dp, mo_coeff, matrix_dSdB_mo, &
141 0 : 0.0_dp, vcd_env%dCB_prime(alpha))
142 :
143 : END DO
144 :
145 : ! The AATs in MFP consist of four terms
146 : !
147 : ! i) (i CB_alpha) * CR_beta * S0
148 : ! ii) (i CB_alpha) * C0 * < mu | d/dR_beta^\lambda nu >
149 : ! iii) -1/c * ε_{alpha gamma delta} * C0 * C0 * < mu | r_delta | d/dR_beta^\lambda nu > * R^mu_gamma
150 : ! iv) -1/c * ε_{alpha gamma delta} * C0 * CR_beta * < mu | r_delta | nu > * R^mu_gamma
151 :
152 : ! iii) -1/c * ε_{alpha gamma delta} * C0 * C0 * < mu | r_delta | d/dR_beta^\lambda nu > * R^mu_gamma
153 : ! +1/c * ε_{alpha gamma delta} * C0 * C0 * < mu | r_delta ∂_beta | nu > * R^mu_gamma * delta_nu^lambda
154 : !
155 : !
156 : ! vcd_env%moments_der_right(delta, beta)%matrix
157 : ! = +/- < mu | r_delta ∂_beta | nu > * (nu == lambda)
158 0 : DO alpha = 1, 3
159 0 : aat_moment_der = 0._dp
160 :
161 0 : DO gamma = 1, 3
162 0 : DO delta = 1, 3
163 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
164 0 : IF (lc_tmp == 0._dp) CYCLE
165 :
166 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
167 0 : vcd_env%moments_der_right(delta, vcd_env%dcdr_env%beta)%matrix)
168 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
169 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
170 0 : gauge_origin=gauge_origin)
171 :
172 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
173 0 : CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
174 :
175 0 : aat_moment_der = aat_moment_der - lc_tmp*tmp_trace
176 : END DO
177 : END DO
178 :
179 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
180 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_moment_der
181 : END DO
182 :
183 : ! And additionally we have the reference dependence
184 : ! + ε_{alpha gamma delta} * P0 * < mu | R^mu_gamma * r_delta ∂_beta | nu > * delta_nu^lambda
185 : ! - ε_{alpha gamma delta} * P0 * < mu | r_delta ∂_beta | nu > * R^G_gamma * delta_nu^lambda
186 : ! - ε_{alpha gamma delta} * P0 * < mu | R^mu_gamma * ∂_beta | nu > * R^G_delta * delta_nu^lambda
187 : ! + ε_{alpha gamma delta} * P0 * < mu | ∂_beta | nu > * R^G_gamma * R^G_delta * delta_nu^lambda
188 0 : DO alpha = 1, 3
189 0 : aat_tmp = 0._dp
190 0 : DO gamma = 1, 3
191 0 : DO delta = 1, 3
192 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
193 0 : IF (lc_tmp == 0._dp) CYCLE
194 :
195 : ! - < mu | r_delta ∂_beta | nu > * R^G_gamma * delta_nu^lambda
196 : ! I have
197 : ! vcd_env%moments_der_right(delta, beta)%matrix = -< mu | r_delta ∂_beta | nu > * (nu == lambda)
198 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%moments_der_right(delta, vcd_env%dcdr_env%beta)%matrix, &
199 0 : mo_coeff, buf, ncol=nmo)
200 0 : CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
201 0 : aat_tmp = aat_tmp + lc_tmp*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
202 :
203 : ! - < mu | R^mu_gamma * ∂_beta | nu > * R^G_delta * delta_nu^lambda
204 : ! I have
205 : ! dipvel_ao(beta) = < a | ∂_beta | b >
206 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dipvel_ao(vcd_env%dcdr_env%beta)%matrix)
207 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
208 0 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE.)
209 : CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
210 0 : sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
211 :
212 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, &
213 0 : ncol=nmo, alpha=1._dp, beta=0._dp)
214 0 : CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
215 0 : aat_tmp = aat_tmp - lc_tmp*tmp_trace*vcd_env%spatial_origin_atom(delta)
216 :
217 : ! + < mu | ∂_beta | nu > * R^G_gamma * R^G_delta * delta_nu^lambda
218 : ! I have
219 : ! dipvel_ao(beta) = < a | ∂_beta | b >
220 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dipvel_ao(vcd_env%dcdr_env%beta)%matrix)
221 : CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
222 0 : sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
223 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, &
224 0 : ncol=nmo, alpha=1._dp, beta=0._dp)
225 0 : CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
226 : aat_tmp = aat_tmp + lc_tmp*tmp_trace &
227 0 : *vcd_env%spatial_origin_atom(delta)*vcd_env%magnetic_origin_atom(gamma)
228 :
229 : END DO
230 : END DO
231 :
232 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
233 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp
234 : END DO
235 :
236 : ! ii) (i CB_alpha) * C0 * < mu | d/dR_beta^\lambda nu >
237 : ! = - (i CB_alpha) * C0 * < mu | ∂_beta | nu > * delta_nu^lambda
238 : !
239 : ! dipvel_ao = + < a | ∂ | b >, so I can take that and multiply with the delta
240 0 : CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(vcd_env%dcdr_env%beta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
241 : CALL hr_mult_by_delta_1d(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, qs_kind_set, "ORB", &
242 0 : sab_all, direction_Or=.TRUE., lambda=vcd_env%dcdr_env%lambda)
243 :
244 0 : DO alpha = 1, 3
245 0 : CALL cp_fm_set_all(buf, 0.0_dp)
246 : aat_linmom = 0.0_dp
247 :
248 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
249 0 : mo_coeff, buf, ncol=nmo, alpha=1._dp, beta=0._dp)
250 0 : CALL cp_fm_trace(buf, vcd_env%dCB(alpha), aat_linmom)
251 :
252 : ! the minus, because dipvel_ao is the derivative wrt to the electronic coordinate,
253 : ! but we need the nuclear coordinate here.
254 0 : aat_linmom = -f_spin*aat_linmom
255 :
256 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
257 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
258 : END DO
259 :
260 0 : DO alpha = 1, 3
261 0 : CALL cp_fm_set_all(buf, 0.0_dp)
262 : aat_linmom = 0.0_dp
263 :
264 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
265 0 : CALL cp_fm_trace(buf, vcd_env%dCB_prime(alpha), aat_linmom)
266 :
267 : ! the minus, because dipvel_ao is the derivative wrt to the electronic coordinate,
268 : ! but we need the nuclear coordinate here.
269 0 : aat_linmom = -f_spin*aat_linmom
270 :
271 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
272 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
273 : END DO
274 :
275 : ! The reference dependent dSdB term is
276 : ! i/2c (- < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma) * R^G_delta
277 0 : DO alpha = 1, 3
278 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
279 :
280 0 : DO gamma = 1, 3
281 0 : DO delta = 1, 3
282 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
283 0 : IF (lc_tmp == 0._dp) CYCLE
284 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
285 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
286 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
287 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
288 :
289 : ! < mu | nu > * R^mu_gamma
290 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
291 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
292 0 : gauge_origin=gauge_origin)
293 :
294 : ! < mu | nu > * R^nu_gamma
295 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
296 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
297 0 : gauge_origin=gauge_origin)
298 :
299 : ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
300 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
301 0 : vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
302 :
303 : ! and accumulate the results
304 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, &
305 : vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
306 0 : 1._dp, vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
307 : END DO
308 : END DO
309 :
310 : ! With the reference dependent overlap derivative in vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix,
311 : ! we now build the corresponding coefficients as buf2 = -0.5 * C0 * S1(ref)
312 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
313 0 : buf, ncol=nmo, alpha=1._dp, beta=0._dp)
314 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
315 : 1.0_dp, mo_coeff, buf, &
316 0 : 0.0_dp, matrix_dSdB_mo)
317 :
318 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
319 : -0.5_dp, mo_coeff, matrix_dSdB_mo, &
320 0 : 0.0_dp, buf2)
321 :
322 : ! vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix still is < a | ∂ | b > * (nu==lambda)
323 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
324 0 : CALL cp_fm_trace(buf, buf2, aat_linmom)
325 :
326 0 : aat_linmom = -f_spin*aat_linmom
327 :
328 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
329 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
330 :
331 : END DO
332 :
333 : ! iv) -1/c * ε_{alpha gamma delta} * C0 * CR_beta * < mu | r_delta | nu > * R^mu_gamma
334 0 : DO alpha = 1, 3
335 0 : aat_moment = 0._dp
336 :
337 0 : DO gamma = 1, 3
338 0 : DO delta = 1, 3
339 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
340 0 : IF (lc_tmp == 0._dp) CYCLE
341 :
342 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
343 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
344 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
345 0 : gauge_origin=gauge_origin)
346 :
347 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
348 0 : CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
349 :
350 0 : aat_moment = aat_moment - lc_tmp*tmp_trace
351 : END DO
352 : END DO
353 :
354 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
355 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_moment
356 : END DO
357 :
358 : ! And additionally we have the reference dependence
359 : ! - ε_{alpha gamma delta} * C0 * CR_beta * < mu | R^mu_gamma r_delta | nu >
360 : ! + ε_{alpha gamma delta} * C0 * CR_beta * < mu | r_delta | nu > * R^G_gamma
361 : ! + ε_{alpha gamma delta} * C0 * CR_beta * < mu | R^mu_gamma | nu > * R^G_delta
362 : ! - ε_{alpha gamma delta} * C0 * CR_beta * < mu | nu > * R^G_gamma * R^G_delta
363 0 : DO alpha = 1, 3
364 0 : aat_moment = 0._dp
365 :
366 0 : DO gamma = 1, 3
367 0 : DO delta = 1, 3
368 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
369 0 : IF (lc_tmp == 0._dp) CYCLE
370 :
371 : ! < mu | r_delta | nu > * R^G_gamma
372 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
373 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
374 0 : CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
375 :
376 0 : aat_moment = aat_moment + lc_tmp*tmp_trace*vcd_env%magnetic_origin_atom(gamma)
377 :
378 : ! < mu | R^mu_gamma | nu > * R^G_delta
379 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
380 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
381 0 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE.)
382 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
383 0 : CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
384 :
385 0 : aat_moment = aat_moment + lc_tmp*tmp_trace*vcd_env%spatial_origin_atom(delta)
386 :
387 : ! < mu | nu > * R^G_gamma * R^G_delta
388 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, mo_coeff, buf, ncol=nmo)
389 0 : CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
390 :
391 : aat_moment = aat_moment + lc_tmp*tmp_trace &
392 0 : *vcd_env%magnetic_origin_atom(gamma)*vcd_env%spatial_origin_atom(delta)
393 : END DO
394 : END DO
395 :
396 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
397 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_moment
398 : END DO
399 :
400 : ! i) 2 * (iCB_alpha) * CR_beta * S0
401 : !
402 0 : DO alpha = 1, 3
403 0 : CALL cp_fm_set_all(buf, 0.0_dp)
404 : aat_overlap = 0.0_dp
405 :
406 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), buf, ncol=nmo)
407 0 : CALL cp_fm_trace(buf, vcd_env%dCB(alpha), aat_overlap)
408 :
409 0 : aat_overlap = f_spin*aat_overlap
410 :
411 : ! dCB_prime * dCR_prime = + iP1
412 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
413 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
414 : END DO
415 :
416 0 : DO alpha = 1, 3
417 0 : CALL cp_fm_set_all(buf, 0.0_dp)
418 : aat_overlap = 0.0_dp
419 :
420 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), buf, ncol=nmo)
421 0 : CALL cp_fm_trace(buf, vcd_env%dCB_prime(alpha), aat_overlap)
422 :
423 0 : aat_overlap = f_spin*aat_overlap
424 :
425 : ! dCB_prime * dCR_prime = + iP1
426 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
427 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
428 : END DO
429 :
430 : ! The reference dependent dSdB term is
431 : ! i/2c (- < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma) * R^G_delta
432 0 : DO alpha = 1, 3
433 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
434 :
435 0 : DO gamma = 1, 3
436 0 : DO delta = 1, 3
437 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
438 0 : IF (lc_tmp == 0._dp) CYCLE
439 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
440 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
441 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
442 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
443 :
444 : ! < mu | nu > * R^mu_gamma
445 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
446 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
447 0 : gauge_origin=gauge_origin)
448 :
449 : ! < mu | nu > * R^nu_gamma
450 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
451 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
452 0 : gauge_origin=gauge_origin)
453 :
454 : ! !! A = alpha*A + beta*B or
455 : ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
456 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
457 0 : vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
458 :
459 : ! and accumulate the results
460 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
461 0 : 1._dp, -vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
462 : END DO
463 : END DO
464 :
465 : ! With the reference dependent overlap derivative in vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix,
466 : ! we now build the corresponding coefficients as buf2 = -0.5 * C0 * S1(ref)
467 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
468 0 : buf, ncol=nmo, alpha=1._dp, beta=0._dp)
469 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
470 : 1.0_dp, mo_coeff, buf, &
471 0 : 0.0_dp, matrix_dSdB_mo)
472 :
473 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
474 : -0.5_dp, mo_coeff, matrix_dSdB_mo, &
475 0 : 0.0_dp, buf2)
476 :
477 0 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%dCR_prime(ispin), buf, ncol=nmo)
478 0 : CALL cp_fm_trace(buf, buf2, aat_overlap)
479 :
480 0 : aat_overlap = f_spin*aat_overlap
481 :
482 : ! dCB_prime * dCR_prime = + iP1
483 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
484 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
485 : END DO
486 :
487 : ! Nuclear contribution
488 0 : CALL get_atomic_kind(particle_set(vcd_env%dcdr_env%lambda)%atomic_kind, kind_number=ikind)
489 0 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
490 0 : IF (.NOT. ghost) THEN
491 0 : DO alpha = 1, 3
492 0 : aat_tmp = 0._dp
493 0 : DO gamma = 1, 3
494 0 : IF (Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta) == 0._dp) CYCLE
495 : aat_tmp = aat_tmp + charge &
496 : *Levi_Civita(alpha, gamma, vcd_env%dcdr_env%beta) &
497 0 : *(particle_set(vcd_env%dcdr_env%lambda)%r(gamma) - vcd_env%magnetic_origin_atom(gamma))
498 :
499 : aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
500 0 : = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_tmp/4.
501 : END DO
502 : END DO
503 : END IF
504 :
505 : END ASSOCIATE
506 :
507 0 : CALL cp_fm_release(buf)
508 0 : CALL cp_fm_release(buf2)
509 0 : CALL cp_fm_release(matrix_dSdB_mo)
510 :
511 0 : CALL timestop(handle)
512 0 : END SUBROUTINE mfp_aat
513 :
514 : ! **************************************************************************************************
515 : !> \brief ...
516 : !> \param vcd_env ...
517 : !> \param qs_env ...
518 : !> \param alpha ...
519 : !> \author Edward Ditler
520 : ! **************************************************************************************************
521 0 : SUBROUTINE mfp_build_operator_gauge_dependent(vcd_env, qs_env, alpha)
522 : TYPE(vcd_env_type) :: vcd_env
523 : TYPE(qs_environment_type), POINTER :: qs_env
524 : INTEGER :: alpha
525 :
526 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mfp_build_operator_gauge_dependent'
527 : INTEGER, PARAMETER :: ispin = 1
528 :
529 : INTEGER :: delta, gamma, handle, nao, nmo
530 : REAL(dp) :: eps_ppnl, lc_tmp
531 : REAL(dp), DIMENSION(3) :: gauge_origin
532 : TYPE(cell_type), POINTER :: cell
533 : TYPE(cp_fm_type) :: buf
534 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
535 : TYPE(dft_control_type), POINTER :: dft_control
536 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
537 0 : POINTER :: sab_all, sap_ppnl
538 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
539 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
540 :
541 : ! --------------------------------------------------------------- !
542 : ! The operator consists of four parts:
543 : ! 1. R^mu x r * H^KS
544 : ! 2. H^KS * R^nu x r
545 : ! 3. [Vnl, (R^mu - Omag) x r]
546 : ! 4. (r - Omag) x ∂
547 : ! and additionally the overlap derivative comes in as
548 : ! 5. -dSdB * C0 * CHC
549 : ! --------------------------------------------------------------- !
550 :
551 0 : CALL timeset(routineN, handle)
552 :
553 : ! Some setup
554 0 : nmo = vcd_env%dcdr_env%nmo(ispin)
555 0 : nao = vcd_env%dcdr_env%nao
556 : ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
557 :
558 : CALL get_qs_env(qs_env=qs_env, &
559 : sab_all=sab_all, &
560 : qs_kind_set=qs_kind_set, &
561 : particle_set=particle_set, &
562 : sap_ppnl=sap_ppnl, &
563 : cell=cell, &
564 : dft_control=dft_control, & ! For the eps_ppnl
565 0 : matrix_ks=matrix_ks)
566 :
567 0 : gauge_origin(:) = vcd_env%magnetic_origin_atom
568 :
569 0 : eps_ppnl = dft_control%qs_control%eps_ppnl
570 :
571 0 : CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
572 0 : CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
573 0 : CALL cp_fm_set_all(buf, 0._dp)
574 :
575 : ! The first two terms in the operator are
576 : ! i/2c sum_{gamma delta} ε_{alpha gamma delta}
577 : ! (R^mu_gamma - O^mag_gamma) * < (r_delta - O^spatial) * H >
578 : ! - (R^nu_gamma - O^mag_gamma) * < H * (r_delta - O^spatial) >)
579 :
580 : ! they expand into
581 : ! i) R^mu_gamma * matrix_rh(delta) - R^nu_gamma * matrix_hr(delta)
582 : ! ii) - O^mag_gamma * (matrix_rh(delta) - matrix_hr(delta))
583 : ! iii) - (R^mu_gamma - R^nu_gamma) * O^spatial_delta * matrix_ks
584 :
585 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
586 0 : DO gamma = 1, 3
587 0 : DO delta = 1, 3
588 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
589 0 : IF (lc_tmp == 0._dp) CYCLE
590 :
591 : ! i)
592 : ! R^mu_gamma * matrix_rh(delta) - R^nu_gamma * matrix_hr(delta)
593 : ! R^mu_gamma * matrix_rh(delta)
594 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
595 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
596 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
597 0 : gauge_origin=gauge_origin)
598 :
599 : ! - R^nu_gamma * matrix_hr(delta)
600 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_hr(ispin, delta)%matrix)
601 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
602 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
603 0 : gauge_origin=gauge_origin)
604 :
605 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
606 0 : vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
607 :
608 : ! and accumulate the results
609 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
610 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
611 :
612 : ! ii)
613 : ! - O^mag_gamma * (matrix_rh(delta) - matrix_hr(delta))
614 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
615 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_hr(ispin, delta)%matrix)
616 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
617 0 : vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
618 0 : CALL dbcsr_scale(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, -vcd_env%magnetic_origin_atom(gamma))
619 :
620 : ! and accumulate the results
621 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
622 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
623 :
624 : ! iii)
625 : ! - (R^mu_gamma - R^nu_gamma) * O^spatial_delta * matrix_ks
626 0 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
627 0 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
628 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
629 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
630 0 : gauge_origin=gauge_origin)
631 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
632 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
633 0 : gauge_origin=gauge_origin)
634 :
635 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
636 0 : 1._dp, -1._dp)
637 0 : CALL dbcsr_scale(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, -vcd_env%spatial_origin_atom(delta))
638 : ! and accumulate the results
639 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
640 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
641 :
642 : END DO
643 : END DO
644 :
645 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
646 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp, beta=1.0_dp)
647 :
648 : ! The third term in the operator is
649 : ! i/2c sum_{gamma delta} ε_{alpha gamma delta}
650 : ! sum_\eta < (R^eta_gamma - O^mag) * [Vnl^eta, r_delta] >
651 : !
652 : ! build_com_vnl_giao calculates
653 : ! matrix_rv(gamma, delta) =
654 : ! sum_\eta < mu | R^eta_gamma * Vnl^eta * r_delta | nu >
655 : ! or the other way around: r_delta * Vnl^eta
656 0 : DO gamma = 1, 3
657 0 : DO delta = 1, 3
658 0 : CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, 0._dp)
659 0 : CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, 0._dp)
660 : END DO
661 : END DO
662 : ! V * r_delta * R^eta_gamma
663 : CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
664 : eps_ppnl=dft_control%qs_control%eps_ppnl, &
665 : particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp_33, &
666 0 : ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.TRUE.)
667 :
668 : ! r_delta * V * R^eta_gamma
669 : CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
670 : eps_ppnl=dft_control%qs_control%eps_ppnl, &
671 : particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp2_33, &
672 0 : ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.FALSE.)
673 :
674 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
675 0 : DO gamma = 1, 3
676 0 : DO delta = 1, 3
677 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
678 0 : IF (lc_tmp == 0._dp) CYCLE
679 :
680 : ! + V * r_delta * R^eta_gamma
681 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, &
682 0 : 1._dp, lc_tmp)
683 :
684 : ! - r_delta * V * R^eta_gamma
685 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, &
686 0 : 1._dp, -lc_tmp)
687 : END DO
688 : END DO
689 :
690 : ! Same factor 1/2 as above for Hr - rH
691 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
692 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp/2._dp, beta=1.0_dp)
693 :
694 : ! And the O^mag dependent term:
695 : ! - sum_\eta < O^mag * [Vnl^eta, r_delta] >
696 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
697 0 : DO gamma = 1, 3
698 0 : DO delta = 1, 3
699 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
700 0 : IF (lc_tmp == 0._dp) CYCLE
701 :
702 : ! vcd_env%hcom(delta) = - < mu | [V, r_delta] | nu >
703 0 : CALL dbcsr_desymmetrize(vcd_env%hcom(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
704 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
705 0 : 1._dp, vcd_env%magnetic_origin_atom(gamma)*lc_tmp/(2._dp))
706 : END DO
707 : END DO
708 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
709 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
710 :
711 : ! The fourth term in the operator is
712 : ! -i/2c sum_{gamma delta} ε_{alpha gamma delta}
713 : ! (r_gamma - O^mag_gamma) * ∂_delta
714 :
715 : ! It's also the P^1 term of the NVPT AAT with the reversed sign
716 : ! So in the first term alpha=+1
717 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
718 0 : DO gamma = 1, 3
719 0 : DO delta = 1, 3
720 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
721 0 : IF (lc_tmp == 0._dp) CYCLE
722 :
723 : ! moments_der(gamma, delta) = - < a | r_gamma \partial_delta | b >
724 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%moments_der(gamma, delta)%matrix, &
725 0 : 1._dp, lc_tmp/(2._dp))
726 :
727 : END DO
728 : END DO
729 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
730 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
731 :
732 : ! And the O^mag dependent term
733 : ! + i/2c sum_{gamma delta} ε_{alpha gamma delta}
734 : ! O^mag_gamma * ∂_delta
735 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
736 0 : DO gamma = 1, 3
737 0 : DO delta = 1, 3
738 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
739 0 : IF (lc_tmp == 0._dp) CYCLE
740 :
741 : ! dipvel_ao(delta) = + < a | ∂_delta | b >
742 0 : CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
743 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
744 0 : 1._dp, vcd_env%magnetic_origin_atom(gamma)*lc_tmp/(2._dp))
745 : END DO
746 : END DO
747 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
748 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
749 :
750 : ! And the overlap derivative
751 : ! i/2c * sum_{gamma delta} ε_{alpha gamma delta}
752 : ! < r_delta > * (R^mu_gamma - R^nu_gamma)
753 : ! The reference independent part
754 0 : CALL dbcsr_set(vcd_env%matrix_dSdB(alpha)%matrix, 0._dp)
755 :
756 0 : DO gamma = 1, 3
757 0 : DO delta = 1, 3
758 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
759 0 : IF (lc_tmp == 0._dp) CYCLE
760 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
761 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 0._dp)
762 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
763 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
764 :
765 : ! moments * R^mu_gamma
766 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
767 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
768 0 : gauge_origin=gauge_origin)
769 :
770 : ! moments * R^nu_gamma
771 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
772 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
773 0 : gauge_origin=gauge_origin)
774 :
775 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
776 0 : vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
777 :
778 : ! and accumulate the results
779 : CALL dbcsr_add(vcd_env%matrix_dSdB(alpha)%matrix, &
780 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
781 : END DO
782 : END DO
783 :
784 : ! the overlap derivative goes in as -S(1,B) * C0 * E0
785 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdB(alpha)%matrix, mo_coeff, &
786 0 : buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
787 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
788 : -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), & ! the negative sign is here
789 0 : 1.0_dp, vcd_env%op_dB(ispin))
790 :
791 : ! And the reference dependent part of the overlap
792 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
793 0 : DO gamma = 1, 3
794 0 : DO delta = 1, 3
795 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
796 0 : IF (lc_tmp == 0._dp) CYCLE
797 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
798 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
799 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
800 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
801 :
802 : ! < mu | nu > * R^mu_gamma
803 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
804 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
805 0 : gauge_origin=gauge_origin)
806 :
807 : ! < mu | nu > * R^nu_gamma
808 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
809 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
810 0 : gauge_origin=gauge_origin)
811 :
812 : ! !! A = alpha*A + beta*B or
813 : ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
814 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
815 0 : vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
816 :
817 : ! and accumulate the results
818 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
819 0 : 1._dp, vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
820 : END DO
821 : END DO
822 :
823 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
824 0 : buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
825 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
826 : -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), & ! again, the negative sign is here
827 0 : 1.0_dp, vcd_env%op_dB(ispin))
828 :
829 : ! CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
830 0 : CALL cp_fm_release(buf)
831 :
832 : END ASSOCIATE
833 :
834 : ! We have built op_dB but plug -op_dB into the linres_solver
835 : ! CALL cp_fm_scale(-1.0_dp, vcd_env%op_dB(ispin))
836 :
837 0 : CALL timestop(handle)
838 0 : END SUBROUTINE mfp_build_operator_gauge_dependent
839 :
840 : ! **************************************************************************************************
841 : !> \brief ...
842 : !> \param vcd_env ...
843 : !> \param qs_env ...
844 : !> \param alpha ...
845 : !> \author Edward Ditler
846 : ! **************************************************************************************************
847 0 : SUBROUTINE mfp_build_operator_gauge_independent(vcd_env, qs_env, alpha)
848 : TYPE(vcd_env_type) :: vcd_env
849 : TYPE(qs_environment_type), POINTER :: qs_env
850 : INTEGER :: alpha
851 :
852 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mfp_build_operator_gauge_independent'
853 : INTEGER, PARAMETER :: ispin = 1
854 :
855 : INTEGER :: delta, gamma, handle, nao, nmo
856 : REAL(dp) :: eps_ppnl, lc_tmp
857 : REAL(dp), DIMENSION(3) :: gauge_origin
858 : TYPE(cell_type), POINTER :: cell
859 : TYPE(cp_fm_type) :: buf
860 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
861 : TYPE(dft_control_type), POINTER :: dft_control
862 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
863 0 : POINTER :: sab_all, sap_ppnl
864 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
865 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
866 :
867 : ! --------------------------------------------------------------- !
868 : ! The operator consists of three parts:
869 : ! 1. (R^mu - R^nu) x r * H_KS
870 : ! 2. [Vnl, (R^mu - R^nu) x r]
871 : ! 3. (r - R^nu) x ∂
872 : ! and additionally the overlap derivative comes in as
873 : ! 4. -dSdB * C0 * CHC
874 : ! --------------------------------------------------------------- !
875 :
876 0 : CALL timeset(routineN, handle)
877 :
878 : ! Some setup
879 0 : nmo = vcd_env%dcdr_env%nmo(ispin)
880 0 : nao = vcd_env%dcdr_env%nao
881 : ASSOCIATE (mo_coeff => vcd_env%dcdr_env%mo_coeff(ispin))
882 :
883 : CALL get_qs_env(qs_env=qs_env, &
884 : sab_all=sab_all, &
885 : qs_kind_set=qs_kind_set, &
886 : particle_set=particle_set, &
887 : sap_ppnl=sap_ppnl, &
888 : cell=cell, &
889 : dft_control=dft_control, & ! For the eps_ppnl
890 0 : matrix_ks=matrix_ks)
891 :
892 0 : gauge_origin(:) = 0._dp
893 :
894 0 : eps_ppnl = dft_control%qs_control%eps_ppnl
895 :
896 0 : CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
897 0 : CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
898 0 : CALL cp_fm_set_all(buf, 0._dp)
899 :
900 : ! The first term in the operator is
901 : ! i/2c sum_{gamma delta} ε_{alpha gamma delta}
902 : ! (R^mu_gamma - R^nu_gamma) < (r_delta - Ospatial) H >
903 :
904 : ! In vcd_second_prepare we build
905 : ! vcd_env%matrix_hr(1, delta) = < mu | H * r_delta | nu >
906 : ! and vcd_env%matrix_rh(1, delta) = < mu | r_delta * H | nu >
907 :
908 : ! So we need
909 : ! (R^mu_gamma - R^nu_gamma) * vcd_env%matrix_rh
910 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
911 0 : DO gamma = 1, 3
912 0 : DO delta = 1, 3
913 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
914 0 : IF (lc_tmp == 0._dp) CYCLE
915 :
916 : ! The reference independent part
917 : ! vcd_env%matrix_rh * R^mu_gamma
918 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
919 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
920 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
921 0 : gauge_origin=gauge_origin)
922 :
923 : ! - vcd_env%matrix_rh * R^nu_gamma
924 0 : CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
925 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
926 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
927 0 : gauge_origin=gauge_origin)
928 :
929 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
930 0 : vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
931 :
932 : ! and accumulate the results
933 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
934 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
935 :
936 : ! And the Ospatial dependent part
937 : ! - (R^mu_gamma - R^nu_gamma) * Ospatial * matrix_ks
938 0 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
939 0 : CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
940 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
941 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
942 0 : gauge_origin=gauge_origin)
943 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
944 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
945 0 : gauge_origin=gauge_origin)
946 :
947 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
948 0 : 1._dp, -1._dp)
949 0 : CALL dbcsr_scale(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, -vcd_env%spatial_origin_atom(delta))
950 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
951 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
952 :
953 : END DO
954 : END DO
955 :
956 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
957 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp, beta=1.0_dp)
958 :
959 : ! The second term in the operator is
960 : ! i/2c sum_{gamma delta} ε_{alpha gamma delta}
961 : ! sum_\eta < (R^eta_gamma - R^nu_gamma) * [Vnl^eta, r_delta] >
962 : !
963 : ! build_com_vnl_giao calculates
964 : ! matrix_rv(gamma, delta) =
965 : ! sum_\eta < mu | R^eta_gamma * Vnl^eta * r_delta | nu >
966 : ! or the other way around: r_delta * Vnl^eta
967 0 : DO gamma = 1, 3
968 0 : DO delta = 1, 3
969 0 : CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, 0._dp)
970 0 : CALL dbcsr_set(vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, 0._dp)
971 : END DO
972 : END DO
973 : ! V * r_delta * R^eta_gamma
974 : CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
975 : eps_ppnl=dft_control%qs_control%eps_ppnl, &
976 : particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp_33, &
977 0 : ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.TRUE.)
978 :
979 : ! r_delta * V * R^eta_gamma
980 : CALL build_com_vnl_giao(qs_kind_set=qs_kind_set, sab_all=sab_all, sap_ppnl=sap_ppnl, &
981 : eps_ppnl=dft_control%qs_control%eps_ppnl, &
982 : particle_set=particle_set, matrix_rv=vcd_env%matrix_nosym_temp2_33, &
983 0 : ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_Or=.FALSE.)
984 :
985 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
986 0 : DO gamma = 1, 3
987 0 : DO delta = 1, 3
988 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
989 0 : IF (lc_tmp == 0._dp) CYCLE
990 :
991 : ! + V * r_delta * R^eta_gamma
992 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
993 : vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, &
994 0 : 1._dp, lc_tmp)
995 :
996 : ! - r_delta * V * R^eta_gamma
997 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, &
998 : vcd_env%matrix_nosym_temp2_33(gamma, delta)%matrix, &
999 0 : 1._dp, -lc_tmp)
1000 : END DO
1001 : END DO
1002 :
1003 : ! Same factor 1/2 as above for Hr - rH
1004 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
1005 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1._dp/2._dp, beta=1.0_dp)
1006 :
1007 : ! The third term in the operator is
1008 : ! -i/2c sum_{gamma delta} ε_{alpha gamma delta}
1009 : ! (r_gamma - R^nu_gamma) * ∂_delta
1010 :
1011 : ! It's also the P^1 term of the NVPT AAT with the reversed sign
1012 : ! So in the first term alpha=+1
1013 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
1014 0 : DO gamma = 1, 3
1015 0 : DO delta = 1, 3
1016 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
1017 0 : IF (lc_tmp == 0._dp) CYCLE
1018 :
1019 : ! moments_der(gamma, delta) = - < a | r_gamma \partial_delta | b >
1020 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%moments_der(gamma, delta)%matrix, &
1021 0 : 1._dp, lc_tmp/(2._dp))
1022 :
1023 : END DO
1024 : END DO
1025 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
1026 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
1027 :
1028 : ! And the R^nu dependent term
1029 : ! + i/2c sum_{gamma delta} ε_{alpha gamma delta}
1030 : ! R^nu_gamma * ∂_delta
1031 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
1032 0 : DO gamma = 1, 3
1033 0 : DO delta = 1, 3
1034 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
1035 0 : IF (lc_tmp == 0._dp) CYCLE
1036 :
1037 : ! dipvel_ao(delta) = + < a | ∂_delta | b >
1038 0 : CALL dbcsr_desymmetrize(vcd_env%dipvel_ao(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
1039 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
1040 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
1041 0 : gauge_origin=gauge_origin)
1042 :
1043 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
1044 0 : 1._dp, lc_tmp/(2._dp))
1045 : END DO
1046 : END DO
1047 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, mo_coeff, &
1048 0 : vcd_env%op_dB(ispin), ncol=nmo, alpha=1.0_dp, beta=1.0_dp)
1049 :
1050 : ! And the overlap derivative
1051 : ! i/2c * sum_{gamma delta} ε_{alpha gamma delta}
1052 : ! < r_delta > * (R^mu_gamma - R^nu_gamma)
1053 : ! The reference independent part
1054 0 : CALL dbcsr_set(vcd_env%matrix_dSdB(alpha)%matrix, 0._dp)
1055 :
1056 0 : DO gamma = 1, 3
1057 0 : DO delta = 1, 3
1058 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
1059 0 : IF (lc_tmp == 0._dp) CYCLE
1060 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
1061 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 0._dp)
1062 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
1063 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix)
1064 :
1065 : ! moments * R^mu_gamma
1066 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
1067 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
1068 0 : gauge_origin=gauge_origin)
1069 :
1070 : ! moments * R^nu_gamma
1071 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, &
1072 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
1073 0 : gauge_origin=gauge_origin)
1074 :
1075 : ! !! A = alpha*A + beta*B or
1076 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
1077 0 : vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
1078 :
1079 : ! and accumulate the results
1080 : CALL dbcsr_add(vcd_env%matrix_dSdB(alpha)%matrix, &
1081 0 : vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 1._dp, lc_tmp/2._dp)
1082 : END DO
1083 : END DO
1084 :
1085 : ! the overlap derivative goes in as -S(1,B) * C0 * E0
1086 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%matrix_dSdB(alpha)%matrix, mo_coeff, &
1087 0 : buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1088 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
1089 : -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), & ! the negative sign is here
1090 0 : 1.0_dp, vcd_env%op_dB(ispin))
1091 :
1092 : ! And the reference dependent part of the overlap
1093 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
1094 0 : DO gamma = 1, 3
1095 0 : DO delta = 1, 3
1096 0 : lc_tmp = Levi_Civita(alpha, gamma, delta)
1097 0 : IF (lc_tmp == 0._dp) CYCLE
1098 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
1099 0 : CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
1100 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
1101 0 : CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix)
1102 :
1103 : ! < mu | nu > * R^mu_gamma
1104 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
1105 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.FALSE., &
1106 0 : gauge_origin=gauge_origin)
1107 :
1108 : ! < mu | nu > * R^nu_gamma
1109 : CALL multiply_by_position(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, &
1110 : qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.TRUE., &
1111 0 : gauge_origin=gauge_origin)
1112 :
1113 : ! !! A = alpha*A + beta*B or
1114 : ! - < mu | nu > * R^mu_gamma + < mu | nu > * R^nu_gamma
1115 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
1116 0 : vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, -1._dp, +1._dp)
1117 :
1118 : ! and accumulate the results
1119 : CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, &
1120 0 : 1._dp, vcd_env%spatial_origin_atom(delta)*lc_tmp/2._dp)
1121 : END DO
1122 : END DO
1123 :
1124 : CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, mo_coeff, &
1125 0 : buf, ncol=nmo, alpha=1.0_dp, beta=0.0_dp)
1126 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
1127 : -1.0_dp, buf, vcd_env%dcdr_env%chc(ispin), & ! again, the negative sign is here
1128 0 : 1.0_dp, vcd_env%op_dB(ispin))
1129 :
1130 : ! CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
1131 0 : CALL cp_fm_release(buf)
1132 :
1133 : END ASSOCIATE
1134 :
1135 : ! We have built op_dB but plug -op_dB into the linres_solver
1136 : ! CALL cp_fm_scale(-1.0_dp, vcd_env%op_dB(ispin))
1137 :
1138 0 : CALL timestop(handle)
1139 0 : END SUBROUTINE mfp_build_operator_gauge_independent
1140 :
1141 : ! *****************************************************************************
1142 : !> \brief Get the dC/dB using the vcd_env%op_dB
1143 : !> \param vcd_env ...
1144 : !> \param p_env ...
1145 : !> \param qs_env ...
1146 : !> \param alpha ...
1147 : !> \author Edward Ditler
1148 : ! **************************************************************************************************
1149 0 : SUBROUTINE mfp_response(vcd_env, p_env, qs_env, alpha)
1150 :
1151 : TYPE(vcd_env_type) :: vcd_env
1152 : TYPE(qs_p_env_type) :: p_env
1153 : TYPE(qs_environment_type), POINTER :: qs_env
1154 : INTEGER, INTENT(IN) :: alpha
1155 :
1156 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mfp_response'
1157 :
1158 : INTEGER :: handle, output_unit
1159 : LOGICAL :: failure, should_stop
1160 0 : TYPE(cp_fm_type), DIMENSION(1) :: h1_psi0, psi1
1161 : TYPE(cp_logger_type), POINTER :: logger
1162 : TYPE(dft_control_type), POINTER :: dft_control
1163 : TYPE(linres_control_type), POINTER :: linres_control
1164 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1165 : TYPE(section_vals_type), POINTER :: lr_section, vcd_section
1166 :
1167 0 : CALL timeset(routineN, handle)
1168 0 : failure = .FALSE.
1169 :
1170 0 : NULLIFY (linres_control, lr_section, logger)
1171 :
1172 : CALL get_qs_env(qs_env=qs_env, &
1173 : dft_control=dft_control, &
1174 : linres_control=linres_control, &
1175 0 : mos=mos)
1176 :
1177 0 : logger => cp_get_default_logger()
1178 0 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
1179 : vcd_section => section_vals_get_subs_vals(qs_env%input, &
1180 0 : "PROPERTIES%LINRES%VCD")
1181 :
1182 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
1183 0 : extension=".linresLog")
1184 0 : IF (output_unit > 0) THEN
1185 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
1186 0 : "*** Self consistent optimization of the magnetic response wavefunction ***"
1187 : END IF
1188 :
1189 : ! allocate the vectors
1190 : ASSOCIATE (psi0_order => vcd_env%dcdr_env%mo_coeff)
1191 0 : CALL cp_fm_create(psi1(1), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
1192 0 : CALL cp_fm_create(h1_psi0(1), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
1193 :
1194 : ! Restart
1195 0 : IF (linres_control%linres_restart) THEN
1196 0 : CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, alpha, "dCdB")
1197 : ELSE
1198 0 : CALL cp_fm_set_all(psi1(1), 0.0_dp)
1199 : END IF
1200 :
1201 0 : IF (output_unit > 0) THEN
1202 : WRITE (output_unit, *) &
1203 0 : "Response to the perturbation operator referring to the magnetic field in "//ACHAR(alpha + 119)
1204 : END IF
1205 :
1206 : ! First response to get dCR
1207 : ! (H0-E0) psi1 = (H1-E1) psi0
1208 : ! psi1 = the perturbed wavefunction
1209 : ! h1_psi0 = (H1-E1)
1210 : ! psi0_order = the unperturbed wavefunction
1211 : ! Second response to get dCB
1212 0 : CALL cp_fm_set_all(vcd_env%dCB(alpha), 0.0_dp)
1213 0 : CALL cp_fm_set_all(h1_psi0(1), 0.0_dp)
1214 0 : CALL cp_fm_to_fm(vcd_env%op_dB(1), h1_psi0(1))
1215 :
1216 0 : linres_control%lr_triplet = .FALSE. ! we do singlet response
1217 0 : linres_control%do_kernel = .FALSE. ! no coupled response since imaginary perturbation
1218 0 : linres_control%converged = .FALSE.
1219 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
1220 0 : output_unit, should_stop)
1221 0 : CALL cp_fm_to_fm(psi1(1), vcd_env%dCB(alpha))
1222 :
1223 : ! Write the new result to the restart file
1224 0 : IF (linres_control%linres_restart) THEN
1225 0 : CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, alpha, "dCdB")
1226 : END IF
1227 :
1228 : ! clean up
1229 0 : CALL cp_fm_release(psi1(1))
1230 0 : CALL cp_fm_release(h1_psi0(1))
1231 : END ASSOCIATE
1232 :
1233 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1234 0 : "PRINT%PROGRAM_RUN_INFO")
1235 :
1236 0 : CALL timestop(handle)
1237 0 : END SUBROUTINE mfp_response
1238 :
1239 : ! **************************************************************************************************
1240 : !> \brief Take matrix < mu | ^O^ | nu > and multiply the blocks with the positions of the basis functions.
1241 : !> The matrix consists of blocks (iatom, jatom) corresponding to (mu, nu).
1242 : !> With basis_function_nu = .TRUE. we have to multiply each block by particle_set(jatom)%r(delta)
1243 : !> With basis_function_nu = .FALSE. we have to multiply each block by particle_set(iatom)%r(delta)
1244 : !> \param matrix The matrix we are operating on
1245 : !> \param qs_kind_set Needed to set up the basis_sets
1246 : !> \param particle_set Needed for the coordinates
1247 : !> \param basis_type Needed to set up the basis_sets
1248 : !> \param sab_nl The supplied neighborlist
1249 : !> \param direction The index of the nuclear position we are multiplying by
1250 : !> \param basis_function_nu Determines whether we multiply by R^nu or R^mu
1251 : !> \param gauge_origin The gauge origin, if we do distributed gauge
1252 : !> \author Edward Ditler
1253 : ! **************************************************************************************************
1254 0 : SUBROUTINE multiply_by_position(matrix, qs_kind_set, particle_set, basis_type, sab_nl, &
1255 : direction, basis_function_nu, gauge_origin)
1256 :
1257 : TYPE(dbcsr_type) :: matrix
1258 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1259 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1260 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
1261 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1262 : POINTER :: sab_nl
1263 : INTEGER :: direction
1264 : LOGICAL :: basis_function_nu
1265 : REAL(dp), DIMENSION(3), OPTIONAL :: gauge_origin
1266 :
1267 : CHARACTER(len=*), PARAMETER :: routineN = 'multiply_by_position'
1268 :
1269 : INTEGER :: handle, iatom, icol, ikind, inode, irow, &
1270 : jatom, jkind, last_jatom, mepos, &
1271 : nkind, nseta, nsetb, nthread
1272 : INTEGER, DIMENSION(3) :: cell
1273 0 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1274 0 : npgfb, nsgfa, nsgfb
1275 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1276 : LOGICAL :: do_symmetric, found
1277 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1278 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_block, rpgfa, rpgfb, scon_a, &
1279 0 : scon_b, zeta, zetb
1280 0 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1281 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1282 : TYPE(neighbor_list_iterator_p_type), &
1283 0 : DIMENSION(:), POINTER :: nl_iterator
1284 :
1285 0 : CALL timeset(routineN, handle)
1286 :
1287 : ! check for symmetry
1288 0 : CPASSERT(SIZE(sab_nl) > 0)
1289 0 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1290 :
1291 : ! prepare basis set and coordinates
1292 0 : nkind = SIZE(qs_kind_set)
1293 0 : ALLOCATE (basis_set_list(nkind))
1294 0 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
1295 :
1296 : nthread = 1
1297 0 : !$ nthread = omp_get_max_threads()
1298 : ! Iterate of neighbor list
1299 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
1300 :
1301 : !$OMP PARALLEL DEFAULT(NONE) &
1302 : !$OMP SHARED (nthread,nl_iterator, do_symmetric) &
1303 : !$OMP SHARED (matrix,basis_set_list) &
1304 : !$OMP SHARED (basis_function_nu, direction, particle_set, gauge_origin) &
1305 : !$OMP PRIVATE (matrix_block,mepos,ikind,jkind,iatom,jatom,cell) &
1306 : !$OMP PRIVATE (basis_set_a,basis_set_b) &
1307 : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a) &
1308 : !$OMP PRIVATE (zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb) &
1309 0 : !$OMP PRIVATE (zetb, scon_a, scon_b, irow, icol, found, inode, last_jatom)
1310 : mepos = 0
1311 : !$ mepos = omp_get_thread_num()
1312 :
1313 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
1314 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
1315 : iatom=iatom, jatom=jatom, inode=inode)
1316 : basis_set_a => basis_set_list(ikind)%gto_basis_set
1317 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
1318 : basis_set_b => basis_set_list(jkind)%gto_basis_set
1319 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
1320 : ! basis ikind
1321 : ! basis ikind
1322 : first_sgfa => basis_set_a%first_sgf
1323 : la_max => basis_set_a%lmax
1324 : la_min => basis_set_a%lmin
1325 : npgfa => basis_set_a%npgf
1326 : nseta = basis_set_a%nset
1327 : nsgfa => basis_set_a%nsgf_set
1328 : rpgfa => basis_set_a%pgf_radius
1329 : set_radius_a => basis_set_a%set_radius
1330 : scon_a => basis_set_a%scon
1331 : zeta => basis_set_a%zet
1332 : ! basis jkind
1333 : first_sgfb => basis_set_b%first_sgf
1334 : lb_max => basis_set_b%lmax
1335 : lb_min => basis_set_b%lmin
1336 : npgfb => basis_set_b%npgf
1337 : nsetb = basis_set_b%nset
1338 : nsgfb => basis_set_b%nsgf_set
1339 : rpgfb => basis_set_b%pgf_radius
1340 : set_radius_b => basis_set_b%set_radius
1341 : scon_b => basis_set_b%scon
1342 : zetb => basis_set_b%zet
1343 :
1344 : IF (inode == 1) last_jatom = 0
1345 :
1346 : ! this guarantees minimum image convention
1347 : ! anything else would not make sense
1348 : IF (jatom == last_jatom) THEN
1349 : CYCLE
1350 : END IF
1351 :
1352 : last_jatom = jatom
1353 :
1354 : IF (do_symmetric) THEN
1355 : IF (iatom <= jatom) THEN
1356 : irow = iatom
1357 : icol = jatom
1358 : ELSE
1359 : irow = jatom
1360 : icol = iatom
1361 : END IF
1362 : ELSE
1363 : irow = iatom
1364 : icol = jatom
1365 : END IF
1366 :
1367 : NULLIFY (matrix_block)
1368 : CALL dbcsr_get_block_p(matrix, irow, icol, matrix_block, found)
1369 : CPASSERT(found)
1370 :
1371 : IF (PRESENT(gauge_origin)) THEN
1372 : IF (basis_function_nu) THEN
1373 : !$OMP CRITICAL(blockadd)
1374 : matrix_block(:, :) = matrix_block(:, :)*(particle_set(jatom)%r(direction) - gauge_origin(direction))
1375 : !$OMP END CRITICAL(blockadd)
1376 : ELSE
1377 : !$OMP CRITICAL(blockadd)
1378 : matrix_block(:, :) = matrix_block(:, :)*(particle_set(iatom)%r(direction) - gauge_origin(direction))
1379 : !$OMP END CRITICAL(blockadd)
1380 : END IF
1381 : ELSE IF (.NOT. PRESENT(gauge_origin)) THEN
1382 : IF (basis_function_nu) THEN
1383 : !$OMP CRITICAL(blockadd)
1384 : matrix_block(:, :) = matrix_block(:, :)*particle_set(jatom)%r(direction)
1385 : !$OMP END CRITICAL(blockadd)
1386 : ELSE
1387 : !$OMP CRITICAL(blockadd)
1388 : matrix_block(:, :) = matrix_block(:, :)*particle_set(iatom)%r(direction)
1389 : !$OMP END CRITICAL(blockadd)
1390 : END IF
1391 : END IF
1392 : END DO
1393 : !$OMP END PARALLEL
1394 0 : CALL neighbor_list_iterator_release(nl_iterator)
1395 :
1396 : ! Release work storage
1397 0 : DEALLOCATE (basis_set_list)
1398 :
1399 0 : CALL timestop(handle)
1400 :
1401 0 : END SUBROUTINE multiply_by_position
1402 :
1403 : END MODULE qs_mfp
|
