(git:374b731)
Loading...
Searching...
No Matches
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!--------------------------------------------------------------------------------------------------!
7MODULE 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/))
81CONTAINS
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
1403END MODULE qs_mfp
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:74
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_linres_types::vcd_env_type
Definition qs_linres_types.F:302
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_p_env_types::qs_p_env_type
Represent a qs system that is perturbed. Can calculate the linear operator and the rhs of the system ...
Definition qs_p_env_types.F:58