(git:3add494)
qs_mfp.F
Go to the documentation of this file.
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_operations, ONLY: cp_dbcsr_sm_fm_multiply
15  USE cp_fm_basic_linalg, ONLY: cp_fm_trace
16  USE cp_fm_types, ONLY: cp_fm_create,&
17  cp_fm_release,&
18  cp_fm_set_all,&
19  cp_fm_to_fm,&
20  cp_fm_type
21  USE cp_log_handling, ONLY: cp_get_default_logger,&
22  cp_logger_type
23  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
24  cp_print_key_unit_nr
25  USE dbcsr_api, ONLY: dbcsr_add,&
26  dbcsr_copy,&
27  dbcsr_desymmetrize,&
28  dbcsr_get_block_p,&
29  dbcsr_p_type,&
30  dbcsr_scale,&
31  dbcsr_set,&
32  dbcsr_type
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  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  POINTER :: sab_all
107  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
109 
110  CALL timeset(routinen, handle)
111 
112  ! Some setup
113  nmo = vcd_env%dcdr_env%nmo(ispin)
114  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  particle_set=particle_set)
122 
123  ! Set up the matrices
124  CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
125  CALL cp_fm_create(buf2, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
126  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  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  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_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_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  DO alpha = 1, 3
159  aat_moment_der = 0._dp
160 
161  DO gamma = 1, 3
162  DO delta = 1, 3
163  lc_tmp = levi_civita(alpha, gamma, delta)
164  IF (lc_tmp == 0._dp) cycle
165 
166  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
167  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  gauge_origin=gauge_origin)
171 
172  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
173  CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
174 
175  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  = 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  DO alpha = 1, 3
189  aat_tmp = 0._dp
190  DO gamma = 1, 3
191  DO delta = 1, 3
192  lc_tmp = levi_civita(alpha, gamma, delta)
193  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  mo_coeff, buf, ncol=nmo)
200  CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
201  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  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  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  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  ncol=nmo, alpha=1._dp, beta=0._dp)
214  CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
215  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  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  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  ncol=nmo, alpha=1._dp, beta=0._dp)
225  CALL cp_fm_trace(buf, mo_coeff, tmp_trace)
226  aat_tmp = aat_tmp + lc_tmp*tmp_trace &
227  *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  = 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  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  sab_all, direction_or=.true., lambda=vcd_env%dcdr_env%lambda)
243 
244  DO alpha = 1, 3
245  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  mo_coeff, buf, ncol=nmo, alpha=1._dp, beta=0._dp)
250  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  aat_linmom = -f_spin*aat_linmom
255 
256  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
257  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_linmom
258  END DO
259 
260  DO alpha = 1, 3
261  CALL cp_fm_set_all(buf, 0.0_dp)
262  aat_linmom = 0.0_dp
263 
264  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
265  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  aat_linmom = -f_spin*aat_linmom
270 
271  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
272  = 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  DO alpha = 1, 3
278  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
279 
280  DO gamma = 1, 3
281  DO delta = 1, 3
282  lc_tmp = levi_civita(alpha, gamma, delta)
283  IF (lc_tmp == 0._dp) cycle
284  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
285  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
286  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
287  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  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  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  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  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  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_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_dp, buf2)
321 
322  ! vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix still is < a | ∂ | b > * (nu==lambda)
323  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
324  CALL cp_fm_trace(buf, buf2, aat_linmom)
325 
326  aat_linmom = -f_spin*aat_linmom
327 
328  aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) &
329  = 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  DO alpha = 1, 3
335  aat_moment = 0._dp
336 
337  DO gamma = 1, 3
338  DO delta = 1, 3
339  lc_tmp = levi_civita(alpha, gamma, delta)
340  IF (lc_tmp == 0._dp) cycle
341 
342  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  gauge_origin=gauge_origin)
346 
347  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
348  CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
349 
350  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  = 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  DO alpha = 1, 3
364  aat_moment = 0._dp
365 
366  DO gamma = 1, 3
367  DO delta = 1, 3
368  lc_tmp = levi_civita(alpha, gamma, delta)
369  IF (lc_tmp == 0._dp) cycle
370 
371  ! < mu | r_delta | nu > * R^G_gamma
372  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
373  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
374  CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
375 
376  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  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  qs_kind_set, particle_set, "ORB", sab_all, gamma, basis_function_nu=.false.)
382  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, mo_coeff, buf, ncol=nmo)
383  CALL cp_fm_trace(buf, vcd_env%dcdr_env%dCR_prime(ispin), tmp_trace)
384 
385  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  CALL cp_dbcsr_sm_fm_multiply(vcd_env%dcdr_env%matrix_s(1)%matrix, mo_coeff, buf, ncol=nmo)
389  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  *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  = 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  DO alpha = 1, 3
403  CALL cp_fm_set_all(buf, 0.0_dp)
404  aat_overlap = 0.0_dp
405 
406  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  CALL cp_fm_trace(buf, vcd_env%dCB(alpha), aat_overlap)
408 
409  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  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
414  END DO
415 
416  DO alpha = 1, 3
417  CALL cp_fm_set_all(buf, 0.0_dp)
418  aat_overlap = 0.0_dp
419 
420  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  CALL cp_fm_trace(buf, vcd_env%dCB_prime(alpha), aat_overlap)
422 
423  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  = 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  DO alpha = 1, 3
433  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
434 
435  DO gamma = 1, 3
436  DO delta = 1, 3
437  lc_tmp = levi_civita(alpha, gamma, delta)
438  IF (lc_tmp == 0._dp) cycle
439  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
440  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
441  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
442  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  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  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  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  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  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_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_dp, buf2)
476 
477  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  CALL cp_fm_trace(buf, buf2, aat_overlap)
479 
480  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  = aat_atom(vcd_env%dcdr_env%beta, alpha, vcd_env%dcdr_env%lambda) + aat_overlap
485  END DO
486 
487  ! Nuclear contribution
488  CALL get_atomic_kind(particle_set(vcd_env%dcdr_env%lambda)%atomic_kind, kind_number=ikind)
489  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
490  IF (.NOT. ghost) THEN
491  DO alpha = 1, 3
492  aat_tmp = 0._dp
493  DO gamma = 1, 3
494  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  *(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  = 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  CALL cp_fm_release(buf)
508  CALL cp_fm_release(buf2)
509  CALL cp_fm_release(matrix_dsdb_mo)
510 
511  CALL timestop(handle)
512  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  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  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  POINTER :: sab_all, sap_ppnl
538  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
539  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  CALL timeset(routinen, handle)
552 
553  ! Some setup
554  nmo = vcd_env%dcdr_env%nmo(ispin)
555  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  matrix_ks=matrix_ks)
566 
567  gauge_origin(:) = vcd_env%magnetic_origin_atom
568 
569  eps_ppnl = dft_control%qs_control%eps_ppnl
570 
571  CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
572  CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
573  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
586  DO gamma = 1, 3
587  DO delta = 1, 3
588  lc_tmp = levi_civita(alpha, gamma, delta)
589  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  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  gauge_origin=gauge_origin)
598 
599  ! - R^nu_gamma * matrix_hr(delta)
600  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  gauge_origin=gauge_origin)
604 
605  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
606  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  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  CALL dbcsr_copy(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, vcd_env%matrix_rh(ispin, delta)%matrix)
615  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  vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 1._dp, -1._dp)
618  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  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  CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
627  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  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  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  1._dp, -1._dp)
637  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  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  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  DO gamma = 1, 3
657  DO delta = 1, 3
658  CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, 0._dp)
659  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  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  ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_or=.false.)
673 
674  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
675  DO gamma = 1, 3
676  DO delta = 1, 3
677  lc_tmp = levi_civita(alpha, gamma, delta)
678  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  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  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  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
697  DO gamma = 1, 3
698  DO delta = 1, 3
699  lc_tmp = levi_civita(alpha, gamma, delta)
700  IF (lc_tmp == 0._dp) cycle
701 
702  ! vcd_env%hcom(delta) = - < mu | [V, r_delta] | nu >
703  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  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  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
718  DO gamma = 1, 3
719  DO delta = 1, 3
720  lc_tmp = levi_civita(alpha, gamma, delta)
721  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  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  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
736  DO gamma = 1, 3
737  DO delta = 1, 3
738  lc_tmp = levi_civita(alpha, gamma, delta)
739  IF (lc_tmp == 0._dp) cycle
740 
741  ! dipvel_ao(delta) = + < a | ∂_delta | b >
742  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  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  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  CALL dbcsr_set(vcd_env%matrix_dSdB(alpha)%matrix, 0._dp)
755 
756  DO gamma = 1, 3
757  DO delta = 1, 3
758  lc_tmp = levi_civita(alpha, gamma, delta)
759  IF (lc_tmp == 0._dp) cycle
760  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
761  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 0._dp)
762  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
763  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  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  gauge_origin=gauge_origin)
774 
775  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
776  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  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  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  1.0_dp, vcd_env%op_dB(ispin))
790 
791  ! And the reference dependent part of the overlap
792  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
793  DO gamma = 1, 3
794  DO delta = 1, 3
795  lc_tmp = levi_civita(alpha, gamma, delta)
796  IF (lc_tmp == 0._dp) cycle
797  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
798  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
799  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
800  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  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  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  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  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  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  1.0_dp, vcd_env%op_dB(ispin))
828 
829  ! CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
830  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  CALL timestop(handle)
838  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  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  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  POINTER :: sab_all, sap_ppnl
864  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
865  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  CALL timeset(routinen, handle)
877 
878  ! Some setup
879  nmo = vcd_env%dcdr_env%nmo(ispin)
880  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  matrix_ks=matrix_ks)
891 
892  gauge_origin(:) = 0._dp
893 
894  eps_ppnl = dft_control%qs_control%eps_ppnl
895 
896  CALL cp_fm_create(buf, vcd_env%dcdr_env%likemos_fm_struct(ispin)%struct)
897  CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
898  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
911  DO gamma = 1, 3
912  DO delta = 1, 3
913  lc_tmp = levi_civita(alpha, gamma, delta)
914  IF (lc_tmp == 0._dp) cycle
915 
916  ! The reference independent part
917  ! vcd_env%matrix_rh * R^mu_gamma
918  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  gauge_origin=gauge_origin)
922 
923  ! - vcd_env%matrix_rh * R^nu_gamma
924  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  gauge_origin=gauge_origin)
928 
929  CALL dbcsr_add(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, &
930  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  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  CALL dbcsr_desymmetrize(matrix_ks(ispin)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
939  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  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  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  1._dp, -1._dp)
949  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  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  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  DO gamma = 1, 3
968  DO delta = 1, 3
969  CALL dbcsr_set(vcd_env%matrix_nosym_temp_33(gamma, delta)%matrix, 0._dp)
970  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  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  ref_point=[0._dp, 0._dp, 0._dp], cell=cell, direction_or=.false.)
984 
985  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
986  DO gamma = 1, 3
987  DO delta = 1, 3
988  lc_tmp = levi_civita(alpha, gamma, delta)
989  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  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  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  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
1014  DO gamma = 1, 3
1015  DO delta = 1, 3
1016  lc_tmp = levi_civita(alpha, gamma, delta)
1017  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  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  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  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(3)%matrix, 0._dp)
1032  DO gamma = 1, 3
1033  DO delta = 1, 3
1034  lc_tmp = levi_civita(alpha, gamma, delta)
1035  IF (lc_tmp == 0._dp) cycle
1036 
1037  ! dipvel_ao(delta) = + < a | ∂_delta | b >
1038  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  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  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  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  CALL dbcsr_set(vcd_env%matrix_dSdB(alpha)%matrix, 0._dp)
1055 
1056  DO gamma = 1, 3
1057  DO delta = 1, 3
1058  lc_tmp = levi_civita(alpha, gamma, delta)
1059  IF (lc_tmp == 0._dp) cycle
1060  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix, 0._dp)
1061  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp(2)%matrix, 0._dp)
1062  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%moments(delta)%matrix, vcd_env%dcdr_env%matrix_nosym_temp(1)%matrix)
1063  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  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  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  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  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  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  1.0_dp, vcd_env%op_dB(ispin))
1091 
1092  ! And the reference dependent part of the overlap
1093  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(3)%matrix, 0._dp)
1094  DO gamma = 1, 3
1095  DO delta = 1, 3
1096  lc_tmp = levi_civita(alpha, gamma, delta)
1097  IF (lc_tmp == 0._dp) cycle
1098  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix, 0._dp)
1099  CALL dbcsr_set(vcd_env%dcdr_env%matrix_nosym_temp2(2)%matrix, 0._dp)
1100  CALL dbcsr_desymmetrize(vcd_env%dcdr_env%matrix_s(1)%matrix, vcd_env%dcdr_env%matrix_nosym_temp2(1)%matrix)
1101  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  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  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  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  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  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  1.0_dp, vcd_env%op_dB(ispin))
1129 
1130  ! CALL cp_fm_set_all(vcd_env%op_dB(ispin), 0._dp)
1131  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  CALL timestop(handle)
1139  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  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  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1165  TYPE(section_vals_type), POINTER :: lr_section, vcd_section
1166 
1167  CALL timeset(routinen, handle)
1168  failure = .false.
1169 
1170  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  mos=mos)
1176 
1177  logger => cp_get_default_logger()
1178  lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
1179  vcd_section => section_vals_get_subs_vals(qs_env%input, &
1180  "PROPERTIES%LINRES%VCD")
1181 
1182  output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
1183  extension=".linresLog")
1184  IF (output_unit > 0) THEN
1185  WRITE (unit=output_unit, fmt="(T10,A,/)") &
1186  "*** 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  CALL cp_fm_create(psi1(1), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
1192  CALL cp_fm_create(h1_psi0(1), vcd_env%dcdr_env%likemos_fm_struct(1)%struct)
1193 
1194  ! Restart
1195  IF (linres_control%linres_restart) THEN
1196  CALL vcd_read_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, alpha, "dCdB")
1197  ELSE
1198  CALL cp_fm_set_all(psi1(1), 0.0_dp)
1199  END IF
1200 
1201  IF (output_unit > 0) THEN
1202  WRITE (output_unit, *) &
1203  "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  CALL cp_fm_set_all(vcd_env%dCB(alpha), 0.0_dp)
1213  CALL cp_fm_set_all(h1_psi0(1), 0.0_dp)
1214  CALL cp_fm_to_fm(vcd_env%op_dB(1), h1_psi0(1))
1215 
1216  linres_control%lr_triplet = .false. ! we do singlet response
1217  linres_control%do_kernel = .false. ! no coupled response since imaginary perturbation
1218  linres_control%converged = .false.
1219  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
1220  output_unit, should_stop)
1221  CALL cp_fm_to_fm(psi1(1), vcd_env%dCB(alpha))
1222 
1223  ! Write the new result to the restart file
1224  IF (linres_control%linres_restart) THEN
1225  CALL vcd_write_restart(qs_env, lr_section, psi1, vcd_env%dcdr_env%lambda, alpha, "dCdB")
1226  END IF
1227 
1228  ! clean up
1229  CALL cp_fm_release(psi1(1))
1230  CALL cp_fm_release(h1_psi0(1))
1231  END associate
1232 
1233  CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1234  "PRINT%PROGRAM_RUN_INFO")
1235 
1236  CALL timestop(handle)
1237  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  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  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1274  npgfb, nsgfa, nsgfb
1275  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1276  LOGICAL :: do_symmetric, found
1277  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1278  REAL(kind=dp), DIMENSION(:, :), POINTER :: matrix_block, rpgfa, rpgfb, scon_a, &
1279  scon_b, zeta, zetb
1280  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  DIMENSION(:), POINTER :: nl_iterator
1284 
1285  CALL timeset(routinen, handle)
1286 
1287  ! check for symmetry
1288  cpassert(SIZE(sab_nl) > 0)
1289  CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
1290 
1291  ! prepare basis set and coordinates
1292  nkind = SIZE(qs_kind_set)
1293  ALLOCATE (basis_set_list(nkind))
1294  CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
1295 
1296  nthread = 1
1297 !$ nthread = omp_get_max_threads()
1298  ! Iterate of neighbor list
1299  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 !$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  CALL neighbor_list_iterator_release(nl_iterator)
1395 
1396  ! Release work storage
1397  DEALLOCATE (basis_set_list)
1398 
1399  CALL timestop(handle)
1400 
1401  END SUBROUTINE multiply_by_position
1402 
1403 END MODULE qs_mfp
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
basis_set_types
Definition: basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: commutator_rpnl.F:14
commutator_rpnl::build_com_vnl_giao
subroutine, public build_com_vnl_giao(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, matrix_rv, ref_point, cell, direction_Or)
Calculate matrix_rv(gamma, delta) = < R^eta_gamma * Vnl * r_delta > for GIAOs.
Definition: commutator_rpnl.F:1778
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition: qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition: qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_linres_methods
localize wavefunctions linear response scf
Definition: qs_linres_methods.F:15
qs_linres_methods::linres_solver
subroutine, public linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
scf loop to optimize the first order wavefunctions (psi1) given a perturbation as an operator applied...
Definition: qs_linres_methods.F:216
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_mfp
Definition: qs_mfp.F:7
qs_mfp::multiply_by_position
subroutine, public multiply_by_position(matrix, qs_kind_set, particle_set, basis_type, sab_nl, direction, basis_function_nu, gauge_origin)
Take matrix < mu | ^O^ | nu > and multiply the blocks with the positions of the basis functions....
Definition: qs_mfp.F:1256
qs_mfp::mfp_build_operator_gauge_independent
subroutine, public mfp_build_operator_gauge_independent(vcd_env, qs_env, alpha)
...
Definition: qs_mfp.F:848
qs_mfp::mfp_response
subroutine, public mfp_response(vcd_env, p_env, qs_env, alpha)
Get the dC/dB using the vcd_envop_dB.
Definition: qs_mfp.F:1150
qs_mfp::mfp_build_operator_gauge_dependent
subroutine, public mfp_build_operator_gauge_dependent(vcd_env, qs_env, alpha)
...
Definition: qs_mfp.F:522
qs_mfp::mfp_aat
subroutine, public mfp_aat(vcd_env, qs_env)
...
Definition: qs_mfp.F:90
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition: qs_neighbor_list_types.F:993
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition: qs_p_env_types.F:14
qs_vcd_ao
Definition: qs_vcd_ao.F:7
qs_vcd_ao::hr_mult_by_delta_1d
subroutine, public hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_Or)
Apply the operator \delta_\mu^\lambda to zero out all elements of the matrix which don't fulfill the ...
Definition: qs_vcd_ao.F:2579
qs_vcd_utils
Definition: qs_vcd_utils.F:8
qs_vcd_utils::vcd_write_restart
subroutine, public vcd_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_write_restart.
Definition: qs_vcd_utils.F:658
qs_vcd_utils::vcd_read_restart
subroutine, public vcd_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_read_restart.
Definition: qs_vcd_utils.F:508