(git:ed6f26b)
Loading...
Searching...
No Matches
qs_dcdr.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the derivatives of the MO coefficients wrt nuclear coordinates
10!> \author Sandra Luber, Edward Ditler
11! **************************************************************************************************
12
13MODULE qs_dcdr
14
15 USE atomic_kind_types, ONLY: get_atomic_kind
16 USE cell_types, ONLY: cell_type,&
17 pbc
18 USE cp_array_utils, ONLY: cp_2d_r_p_type
19 USE cp_dbcsr_api, ONLY: dbcsr_add,&
20 dbcsr_copy,&
21 dbcsr_desymmetrize,&
22 dbcsr_p_type,&
23 dbcsr_set
24 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
25 dbcsr_allocate_matrix_set,&
26 dbcsr_deallocate_matrix_set
27 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
28 cp_fm_scale_and_add,&
29 cp_fm_trace
30 USE cp_fm_types, ONLY: cp_fm_create,&
31 cp_fm_get_diag,&
32 cp_fm_release,&
33 cp_fm_set_all,&
34 cp_fm_to_fm,&
35 cp_fm_type
36 USE cp_log_handling, ONLY: cp_get_default_logger,&
37 cp_logger_type
38 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
39 cp_print_key_unit_nr
40 USE input_section_types, ONLY: section_vals_get_subs_vals,&
41 section_vals_type
42 USE kinds, ONLY: dp
43 USE molecule_types, ONLY: molecule_of_atom,&
44 molecule_type
45 USE parallel_gemm_api, ONLY: parallel_gemm
46 USE particle_types, ONLY: particle_type
47 USE qs_dcdr_ao, ONLY: apply_op_constant_term,&
48 core_dr,&
49 d_core_charge_density_dr,&
50 d_vhxc_dr,&
51 hr_mult_by_delta_1d,&
52 vhxc_r_perturbed_basis_functions
53 USE qs_dcdr_utils, ONLY: dcdr_read_restart,&
54 dcdr_write_restart,&
55 multiply_localization,&
56 shift_wannier_into_cell
57 USE qs_environment_types, ONLY: get_qs_env,&
58 qs_environment_type
59 USE qs_kind_types, ONLY: get_qs_kind,&
60 qs_kind_type
61 USE qs_linres_methods, ONLY: linres_solver
62 USE qs_linres_types, ONLY: dcdr_env_type,&
63 get_polar_env,&
64 linres_control_type,&
65 polar_env_type
66 USE qs_mo_types, ONLY: get_mo_set,&
67 mo_set_type
68 USE qs_moments, ONLY: build_local_moment_matrix,&
69 dipole_deriv_ao
70 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
71 USE qs_p_env_types, ONLY: qs_p_env_type
72#include "./base/base_uses.f90"
73
74 IMPLICIT NONE
75
76 PRIVATE
77 PUBLIC :: prepare_per_atom, dcdr_response_dr, dcdr_build_op_dr, apt_dr, apt_dr_localization
78
79 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr'
80
81CONTAINS
82
83! **************************************************************************************************
84!> \brief Prepare the environment for a choice of lambda
85!> \param dcdr_env ...
86!> \param qs_env ...
87!> \author Edward Ditler
88! **************************************************************************************************
89 SUBROUTINE prepare_per_atom(dcdr_env, qs_env)
90 TYPE(dcdr_env_type) :: dcdr_env
91 TYPE(qs_environment_type), POINTER :: qs_env
92
93 CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_per_atom'
94
95 INTEGER :: handle, i, ispin, j, natom
96 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
97 POINTER :: sab_all
98 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
99 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
100
101 CALL timeset(routinen, handle)
102
103 NULLIFY (sab_all, qs_kind_set, particle_set)
104 CALL get_qs_env(qs_env=qs_env, &
105 sab_all=sab_all, &
106 qs_kind_set=qs_kind_set, &
107 particle_set=particle_set)
108
109 natom = SIZE(particle_set)
110 IF (dcdr_env%distributed_origin) dcdr_env%ref_point(:) = particle_set(dcdr_env%lambda)%r(:)
111
112 dcdr_env%delta_basis_function = 0._dp
113 dcdr_env%delta_basis_function(:, dcdr_env%lambda) = 1._dp
114
115 ! S matrix
116 ! S1 = - < da/dr | b > * delta_a - < a | db/dr > * delta_b
117
118 ! matrix_s(2:4) are anti-symmetric matrices and contain derivatives wrt. to < a |
119 ! = < da/dR | b > = - < da/dr | b > = < a | db/dr >
120 ! matrix_s1(2:4) = d/dR < a | b >
121 ! and it's built as
122 ! = - matrix_s * delta_b + matrix_s * delta_a
123 ! = - < da/dR | b > * delta_b + < da/dR | b > * delta_a
124 ! = + < da/dr | b > * delta_b - < da/dr | b > * delta_a
125 ! = - < a | db/dr > * delta_b - < da/dr | b > * delta_a
126
127 DO i = 1, 3
128 ! S matrix
129 CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
130 CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_s1(1 + i)%matrix)
131 CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
132
133 CALL hr_mult_by_delta_1d(dcdr_env%matrix_s1(1 + i)%matrix, qs_kind_set, "ORB", &
134 sab_all, dcdr_env%lambda, direction_or=.true.)
135 CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
136 sab_all, dcdr_env%lambda, direction_or=.false.)
137
138 CALL dbcsr_add(dcdr_env%matrix_s1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
139 CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
140
141 ! T matrix
142 CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
143 CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_t1(1 + i)%matrix)
144 CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
145
146 CALL hr_mult_by_delta_1d(dcdr_env%matrix_t1(1 + i)%matrix, qs_kind_set, "ORB", &
147 sab_all, dcdr_env%lambda, direction_or=.true.)
148 CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
149 sab_all, dcdr_env%lambda, direction_or=.false.)
150
151 CALL dbcsr_add(dcdr_env%matrix_t1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
152 CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
153 END DO
154
155 ! Operator:
156 DO ispin = 1, dcdr_env%nspins
157 DO i = 1, 3
158 CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
159 CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
160 CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i)%matrix, 0.0_dp)
161 CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i + 3)%matrix, 0.0_dp)
162 CALL dbcsr_set(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, 0.0_dp)
163 CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
164 END DO
165 END DO
166
167 CALL core_dr(qs_env, dcdr_env) ! dcdr_env%matrix_ppnl_1, hc
168 CALL d_vhxc_dr(qs_env, dcdr_env) ! dcdr_env%matrix_d_vhxc_dR
169 CALL d_core_charge_density_dr(qs_env, dcdr_env) ! dcdr_env%matrix_core_charge_1
170 CALL vhxc_r_perturbed_basis_functions(qs_env, dcdr_env) ! dcdr_env%matrix_vhxc_perturbed_basis
171
172 ! APT:
173 DO i = 1, 3
174 DO j = 1, 3
175 CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0._dp)
176 END DO
177 END DO
178
179 CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, 1, dcdr_env%ref_point)
180
181 CALL timestop(handle)
182 END SUBROUTINE prepare_per_atom
183
184! **************************************************************************************************
185!> \brief Build the operator for the position perturbation
186!> \param dcdr_env ...
187!> \param qs_env ...
188!> \authors Sandra Luber
189!> Edward Ditler
190!> Ravi Kumar
191!> Rangsiman Ketkaew
192! **************************************************************************************************
193 SUBROUTINE dcdr_build_op_dr(dcdr_env, qs_env)
194
195 TYPE(dcdr_env_type) :: dcdr_env
196 TYPE(qs_environment_type), POINTER :: qs_env
197
198 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_build_op_dR'
199 REAL(kind=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
200
201 INTEGER :: handle, ispin, nao, nmo
202 TYPE(cp_fm_type) :: buf
203 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: opdr_sym
204
205 CALL timeset(routinen, handle)
206
207 nao = dcdr_env%nao
208
209 ! allocate matrix for the sum of the perturbation terms of the operator (dbcsr matrix)
210 NULLIFY (opdr_sym)
211 CALL dbcsr_allocate_matrix_set(opdr_sym, 1)
212 ALLOCATE (opdr_sym(1)%matrix)
213 CALL dbcsr_copy(opdr_sym(1)%matrix, dcdr_env%matrix_s1(1)%matrix) ! symmetric
214 CALL dbcsr_set(opdr_sym(1)%matrix, 0.0_dp)
215
216 DO ispin = 1, dcdr_env%nspins
217 nmo = dcdr_env%nmo(ispin)
218
219 CALL apply_op_constant_term(qs_env, dcdr_env) ! dcdr_env%matrix_apply_op_constant
220 ! Hartree and Exchange-Correlation contributions
221 CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_core_charge_1(dcdr_env%beta)%matrix, zero, one)
222 CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_d_vhxc_dR(dcdr_env%beta, ispin)%matrix, one, one)
223 CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_vhxc_perturbed_basis(ispin, dcdr_env%beta)%matrix, one, one)
224
225 ! Core Hamiltonian contributions
226 CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_hc(dcdr_env%beta)%matrix, one, one)
227 CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_ppnl_1(dcdr_env%beta)%matrix, one, one)
228 CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_apply_op_constant(ispin)%matrix, one, one)
229
230 CALL dbcsr_desymmetrize(opdr_sym(1)%matrix, dcdr_env%hamiltonian1(1)%matrix)
231 CALL dbcsr_add(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%matrix_t1(dcdr_env%beta + 1)%matrix, one, one)
232
233 CALL cp_dbcsr_sm_fm_multiply(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%mo_coeff(ispin), &
234 dcdr_env%op_dR(ispin), ncol=nmo)
235
236 ! The overlap derivative terms for the Sternheimer equation
237 ! buf = mo * (-mo * matrix_ks * mo)
238 CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
239 CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
240 -1.0_dp, dcdr_env%mo_coeff(ispin), dcdr_env%chc(ispin), &
241 0.0_dp, buf)
242
243 CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, buf, dcdr_env%op_dR(ispin), &
244 nmo, alpha=1.0_dp, beta=1.0_dp)
245 CALL cp_fm_release(buf)
246
247 ! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
248 CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(ispin))
249
250 IF (dcdr_env%z_matrix_method) THEN
251 CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin))
252 END IF
253
254 END DO
255
256 CALL dbcsr_deallocate_matrix_set(opdr_sym)
257
258 CALL timestop(handle)
259 END SUBROUTINE dcdr_build_op_dr
260
261! **************************************************************************************************
262!> \brief Get the dC/dR by solving the Sternheimer equation, using the op_dR matrix
263!> \param dcdr_env ...
264!> \param p_env ...
265!> \param qs_env ...
266!> \authors SL, ED
267! **************************************************************************************************
268 SUBROUTINE dcdr_response_dr(dcdr_env, p_env, qs_env)
269
270 TYPE(dcdr_env_type) :: dcdr_env
271 TYPE(qs_p_env_type) :: p_env
272 TYPE(qs_environment_type), POINTER :: qs_env
273
274 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_response_dR'
275
276 INTEGER :: handle, ispin, output_unit
277 LOGICAL :: should_stop
278 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
279 TYPE(cp_fm_type), POINTER :: mo_coeff
280 TYPE(cp_logger_type), POINTER :: logger
281 TYPE(linres_control_type), POINTER :: linres_control
282 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
283 TYPE(section_vals_type), POINTER :: lr_section
284
285 CALL timeset(routinen, handle)
286 NULLIFY (linres_control, lr_section, logger)
287
288 CALL get_qs_env(qs_env=qs_env, &
289 linres_control=linres_control, &
290 mos=mos)
291
292 logger => cp_get_default_logger()
293 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
294
295 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
296 extension=".linresLog")
297 IF (output_unit > 0) THEN
298 WRITE (unit=output_unit, fmt="(T10,A,/)") &
299 "*** Self consistent optimization of the response wavefunction ***"
300 END IF
301
302 ! allocate the vectors
303 ALLOCATE (psi0_order(dcdr_env%nspins))
304 ALLOCATE (psi1(dcdr_env%nspins))
305 ALLOCATE (h1_psi0(dcdr_env%nspins))
306
307 DO ispin = 1, dcdr_env%nspins
308 CALL cp_fm_create(psi1(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
309 CALL cp_fm_create(h1_psi0(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
310 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
311 psi0_order(ispin) = mo_coeff
312 END DO
313
314 ! Restart
315 IF (linres_control%linres_restart) THEN
316 CALL dcdr_read_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
317 ELSE
318 DO ispin = 1, dcdr_env%nspins
319 CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
320 END DO
321 END IF
322
323 IF (output_unit > 0) THEN
324 WRITE (output_unit, "(T10,A,I4,A)") &
325 "Response to the perturbation operator referring to atom ", dcdr_env%lambda, &
326 " displaced in "//achar(dcdr_env%beta + 119)
327 END IF
328 DO ispin = 1, dcdr_env%nspins
329 CALL cp_fm_set_all(dcdr_env%dCR(ispin), 0.0_dp)
330 CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), h1_psi0(ispin))
331 END DO
332
333 linres_control%lr_triplet = .false. ! we do singlet response
334 linres_control%do_kernel = .true.
335 linres_control%converged = .false.
336
337 ! Position perturbation to get dCR
338 ! (H0-E0) psi1 = (H1-E1) psi0
339 ! psi1 = the perturbed wavefunction
340 ! h1_psi0 = (H1-E1-S1*\varepsilon)
341 ! psi0_order = the unperturbed wavefunction
342 CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
343 output_unit, should_stop)
344 DO ispin = 1, dcdr_env%nspins
345 CALL cp_fm_to_fm(psi1(ispin), dcdr_env%dCR(ispin))
346 END DO
347
348 ! Write the new result to the restart file
349 IF (linres_control%linres_restart) THEN
350 CALL dcdr_write_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
351 END IF
352
353 ! clean up
354 DO ispin = 1, dcdr_env%nspins
355 CALL cp_fm_release(psi1(ispin))
356 CALL cp_fm_release(h1_psi0(ispin))
357 END DO
358 DEALLOCATE (psi1, h1_psi0, psi0_order)
359 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
360 "PRINT%PROGRAM_RUN_INFO")
361
362 CALL timestop(handle)
363
364 END SUBROUTINE dcdr_response_dr
365
366! **************************************************************************************************
367!> \brief Calculate atomic polar tensor
368!> \param qs_env ...
369!> \param dcdr_env ...
370!> \authors Sandra Luber
371!> Edward Ditler
372!> Ravi Kumar
373!> Rangsiman Ketkaew
374! **************************************************************************************************
375 SUBROUTINE apt_dr(qs_env, dcdr_env)
376 TYPE(qs_environment_type), POINTER :: qs_env
377 TYPE(dcdr_env_type) :: dcdr_env
378
379 CHARACTER(LEN=*), PARAMETER :: routinen = 'apt_dR'
380
381 INTEGER :: alpha, handle, ikind, ispin, nao, nmo
382 LOGICAL :: ghost
383 REAL(dp) :: apt_basis_derivative, &
384 apt_coeff_derivative, charge, f_spin, &
385 temp1, temp2
386 REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
387 TYPE(cp_fm_type) :: overlap1_mo, tmp_fm_like_mos
388 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0, psi1_dberry
389 TYPE(cp_fm_type), POINTER :: mo_coeff
390 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
391 TYPE(polar_env_type), POINTER :: polar_env
392 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
393
394 apt_basis_derivative = 0._dp
395 apt_coeff_derivative = 0._dp
396
397 CALL timeset(routinen, handle)
398
399 NULLIFY (qs_kind_set, particle_set)
400 CALL get_qs_env(qs_env=qs_env, &
401 qs_kind_set=qs_kind_set, &
402 particle_set=particle_set)
403
404 nao = dcdr_env%nao
405 apt_el => dcdr_env%apt_el_dcdr
406 apt_nuc => dcdr_env%apt_nuc_dcdr
407
408 f_spin = 2._dp/dcdr_env%nspins
409
410 DO ispin = 1, dcdr_env%nspins
411 ! Compute S^(1,R)_(ij)
412 CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
413 CALL cp_fm_create(overlap1_mo, dcdr_env%momo_fm_struct(ispin)%struct)
414 nmo = dcdr_env%nmo(ispin)
415 mo_coeff => dcdr_env%mo_coeff(ispin)
416 CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
417 CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
418 CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
419 tmp_fm_like_mos, ncol=nmo)
420 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
421 1.0_dp, mo_coeff, tmp_fm_like_mos, &
422 0.0_dp, overlap1_mo)
423
424 ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
425 ! We get the negative of the coefficients out of the linres solver
426 ! And apply the constant correction due to the overlap derivative.
427 CALL parallel_gemm("N", "N", nao, nmo, nmo, &
428 -0.5_dp, mo_coeff, overlap1_mo, &
429 -1.0_dp, dcdr_env%dCR_prime(ispin))
430 CALL cp_fm_release(overlap1_mo)
431
432 DO alpha = 1, 3
433 IF (.NOT. dcdr_env%z_matrix_method) THEN
434
435 ! FIRST CONTRIBUTION: dCR * moments * mo
436 CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
437 CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
438 CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
439 CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
440 -dcdr_env%ref_point(alpha), 1._dp)
441
442 CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(ispin), &
443 tmp_fm_like_mos, ncol=nmo)
444 CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
445
446 apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
447 apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
448 = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
449 ELSE
450 CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
451 CALL get_polar_env(polar_env=polar_env, psi1_dberry=psi1_dberry, &
452 dberry_psi0=dberry_psi0)
453
454 ! Note that here dcdr_env%dCR_prime contains only occ-occ block contribution,
455 ! dcdr_env%dCR(ispin) is zero because we didn't run response calculation for dcdR.
456
457 CALL cp_fm_trace(dberry_psi0(alpha, ispin), &
458 dcdr_env%dCR_prime(ispin), &
459 temp1)
460
461 CALL cp_fm_trace(dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
462 psi1_dberry(alpha, ispin), &
463 temp2)
464
465 apt_coeff_derivative = temp1 - temp2
466
467 ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
468 ! - apt_coeff_derivative , here the trace is negative to compensate the
469 ! -ve sign in APTs= - 2 Z. M_alpha
470
471 apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
472 apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
473 = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
474 END IF
475
476 ! SECOND CONTRIBUTION: We assemble all combinations of r_i, d(chi)/d(idir)
477 ! difdip contains derivatives with respect to atom dcdr_env%lambda
478 ! difdip(alpha, beta): < a | r_alpha | db/dR_beta >
479 ! Multiply by the MO coefficients
480 CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
481 CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, mo_coeff, &
482 tmp_fm_like_mos, ncol=nmo)
483 CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_basis_derivative)
484
485 ! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
486 apt_basis_derivative = -f_spin*apt_basis_derivative
487 apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
488 apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
489
490 END DO ! alpha
491
492 CALL cp_fm_release(tmp_fm_like_mos)
493 END DO !ispin
494
495 ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
496 CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
497 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
498 IF (.NOT. ghost) THEN
499 apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
500 apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
501 END IF
502
503 ! And deallocate all the things!
504 CALL cp_fm_release(tmp_fm_like_mos)
505 CALL cp_fm_release(overlap1_mo)
506
507 CALL timestop(handle)
508 END SUBROUTINE apt_dr
509
510! **************************************************************************************************
511!> \brief Calculate atomic polar tensor using the localized dipole operator
512!> \param qs_env ...
513!> \param dcdr_env ...
514!> \authors Edward Ditler
515!> Ravi Kumar
516!> Rangsiman Ketkaew
517! **************************************************************************************************
518 SUBROUTINE apt_dr_localization(qs_env, dcdr_env)
519 TYPE(qs_environment_type), POINTER :: qs_env
520 TYPE(dcdr_env_type) :: dcdr_env
521
522 CHARACTER(LEN=*), PARAMETER :: routinen = 'apt_dR_localization'
523
524 INTEGER :: alpha, handle, i, icenter, ikind, ispin, &
525 map_atom, map_molecule, &
526 max_nbr_center, nao, natom, nmo, &
527 nsubset
528 INTEGER, ALLOCATABLE, DIMENSION(:) :: mapping_atom_molecule
529 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mapping_wannier_atom
530 LOGICAL :: ghost
531 REAL(dp) :: apt_basis_derivative, &
532 apt_coeff_derivative, charge, f_spin, &
533 smallest_r, this_factor, tmp_aptcontr, &
534 tmp_r
535 REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements, diagonal_elements2
536 REAL(dp), DIMENSION(3) :: distance, r_shifted
537 REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
538 REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
539 TYPE(cell_type), POINTER :: cell
540 TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
541 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0, psi1_dberry
542 TYPE(cp_fm_type), POINTER :: mo_coeff, overlap1_mo, tmp_fm, &
543 tmp_fm_like_mos, tmp_fm_momo, &
544 tmp_fm_momo2
545 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
546 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
547 TYPE(polar_env_type), POINTER :: polar_env
548 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
549
550 CALL timeset(routinen, handle)
551
552 NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
553
554 CALL get_qs_env(qs_env=qs_env, &
555 qs_kind_set=qs_kind_set, &
556 particle_set=particle_set, &
557 molecule_set=molecule_set, &
558 cell=cell)
559
560 nsubset = SIZE(molecule_set)
561 natom = SIZE(particle_set)
562 apt_el => dcdr_env%apt_el_dcdr
563 apt_nuc => dcdr_env%apt_nuc_dcdr
564 apt_subset => dcdr_env%apt_el_dcdr_per_subset
565 apt_center => dcdr_env%apt_el_dcdr_per_center
566
567 ! Map wannier functions to atoms
568 IF (dcdr_env%nspins == 1) THEN
569 max_nbr_center = dcdr_env%nbr_center(1)
570 ELSE
571 max_nbr_center = max(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
572 END IF
573 ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
574 ALLOCATE (mapping_atom_molecule(natom))
575 centers_set => dcdr_env%centers_set
576
577 DO ispin = 1, dcdr_env%nspins
578 DO icenter = 1, dcdr_env%nbr_center(ispin)
579 ! For every center we check which atom is closest
580 CALL shift_wannier_into_cell(r=centers_set(ispin)%array(1:3, icenter), &
581 cell=cell, &
582 r_shifted=r_shifted)
583
584 smallest_r = huge(0._dp)
585 DO i = 1, natom
586 distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
587 tmp_r = sum(distance**2)
588 IF (tmp_r < smallest_r) THEN
589 mapping_wannier_atom(icenter, ispin) = i
590 smallest_r = tmp_r
591 END IF
592 END DO
593 END DO
594
595 ! Map atoms to molecules
596 CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
597 IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
598 DO icenter = 1, dcdr_env%nbr_center(ispin)
599 map_atom = mapping_wannier_atom(icenter, ispin)
600 map_molecule = mapping_atom_molecule(map_atom)
601 END DO
602 END IF
603 END DO !ispin
604
605 nao = dcdr_env%nao
606 f_spin = 2._dp/dcdr_env%nspins
607
608 DO ispin = 1, dcdr_env%nspins
609 ! Compute S^(1,R)_(ij)
610
611 ALLOCATE (tmp_fm_like_mos)
612 ALLOCATE (overlap1_mo)
613 CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
614 CALL cp_fm_create(overlap1_mo, dcdr_env%momo_fm_struct(ispin)%struct)
615 nmo = dcdr_env%nmo(ispin)
616 mo_coeff => dcdr_env%mo_coeff(ispin)
617 CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
618 CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
619 CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
620 tmp_fm_like_mos, ncol=nmo)
621 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
622 1.0_dp, mo_coeff, tmp_fm_like_mos, &
623 0.0_dp, overlap1_mo)
624
625 ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
626 ! We get the negative of the coefficients out of the linres solver
627 ! And apply the constant correction due to the overlap derivative.
628 CALL parallel_gemm("N", "N", nao, nmo, nmo, &
629 -0.5_dp, mo_coeff, overlap1_mo, &
630 -1.0_dp, dcdr_env%dCR_prime(ispin))
631 CALL cp_fm_release(overlap1_mo)
632
633 ALLOCATE (diagonal_elements(nmo))
634 ALLOCATE (diagonal_elements2(nmo))
635
636 ! Allocate temporary matrices
637 ALLOCATE (tmp_fm)
638 ALLOCATE (tmp_fm_momo)
639 ALLOCATE (tmp_fm_momo2)
640 CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
641 CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
642 CALL cp_fm_create(tmp_fm_momo2, dcdr_env%momo_fm_struct(ispin)%struct)
643
644 ! FIRST CONTRIBUTION: dCR * moments * mo
645 this_factor = -2._dp*f_spin
646 DO alpha = 1, 3
647 IF (.NOT. dcdr_env%z_matrix_method) THEN
648
649 DO icenter = 1, dcdr_env%nbr_center(ispin)
650 CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
651 CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
652 ref_point=centers_set(ispin)%array(1:3, icenter))
653 CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
654 mo_coeff=dcdr_env%dCR_prime(ispin), work=tmp_fm, nmo=nmo, &
655 icenter=icenter, &
656 res=tmp_fm_like_mos)
657 END DO
658
659 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
660 1.0_dp, mo_coeff, tmp_fm_like_mos, &
661 0.0_dp, tmp_fm_momo)
662 CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
663
664 ELSE
665 CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
666 CALL get_polar_env(polar_env=polar_env, psi1_dberry=psi1_dberry, &
667 dberry_psi0=dberry_psi0)
668
669 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
670 1.0_dp, dcdr_env%dCR_prime(ispin), dberry_psi0(alpha, ispin), &
671 0.0_dp, tmp_fm_momo)
672 CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
673
674 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
675 1.0_dp, dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
676 psi1_dberry(alpha, ispin), 0.0_dp, tmp_fm_momo2)
677 CALL cp_fm_get_diag(tmp_fm_momo2, diagonal_elements2)
678
679 diagonal_elements(:) = diagonal_elements(:) - diagonal_elements2(:)
680 END IF
681
682 DO icenter = 1, dcdr_env%nbr_center(ispin)
683 map_atom = mapping_wannier_atom(icenter, ispin)
684 map_molecule = mapping_atom_molecule(map_atom)
685 tmp_aptcontr = this_factor*diagonal_elements(icenter)
686
687 apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
688 = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
689
690 apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
691 = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
692 END DO
693
694 apt_coeff_derivative = this_factor*sum(diagonal_elements)
695 apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
696 = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
697 END DO
698
699 ! SECOND CONTRIBUTION: We assemble all combinations of r_i, dphi/d(idir)
700 ! build part with AOs differentiated with respect to nuclear coordinates
701 ! difdip contains derivatives with respect to atom dcdr_env%lambda
702 ! difdip(alpha, beta): < a | r_alpha | d b/dR_beta >
703 this_factor = -f_spin
704 DO alpha = 1, 3
705 DO icenter = 1, dcdr_env%nbr_center(ispin)
706 ! Build the AO matrix with the right wannier center as reference point
707 CALL dbcsr_set(dcdr_env%matrix_difdip(1, dcdr_env%beta)%matrix, 0._dp)
708 CALL dbcsr_set(dcdr_env%matrix_difdip(2, dcdr_env%beta)%matrix, 0._dp)
709 CALL dbcsr_set(dcdr_env%matrix_difdip(3, dcdr_env%beta)%matrix, 0._dp)
710 CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, &
711 1, centers_set(ispin)%array(1:3, icenter))
712 CALL multiply_localization(ao_matrix=dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, &
713 mo_coeff=mo_coeff, work=tmp_fm, nmo=nmo, &
714 icenter=icenter, &
715 res=tmp_fm_like_mos)
716 END DO ! icenter
717
718 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
719 1.0_dp, mo_coeff, tmp_fm_like_mos, &
720 0.0_dp, tmp_fm_momo)
721 CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
722
723 DO icenter = 1, dcdr_env%nbr_center(ispin)
724 map_atom = mapping_wannier_atom(icenter, ispin)
725 map_molecule = mapping_atom_molecule(map_atom)
726 tmp_aptcontr = this_factor*diagonal_elements(icenter)
727
728 apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
729 = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
730
731 apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
732 = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
733 END DO
734
735 ! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
736 apt_basis_derivative = this_factor*sum(diagonal_elements)
737
738 apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
739 = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
740
741 END DO ! alpha
742 DEALLOCATE (diagonal_elements)
743 DEALLOCATE (diagonal_elements2)
744
745 CALL cp_fm_release(tmp_fm)
746 CALL cp_fm_release(tmp_fm_like_mos)
747 CALL cp_fm_release(tmp_fm_momo)
748 CALL cp_fm_release(tmp_fm_momo2)
749 DEALLOCATE (overlap1_mo)
750 DEALLOCATE (tmp_fm)
751 DEALLOCATE (tmp_fm_like_mos)
752 DEALLOCATE (tmp_fm_momo)
753 DEALLOCATE (tmp_fm_momo2)
754 END DO !ispin
755
756 ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
757 CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
758 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
759 IF (.NOT. ghost) THEN
760 apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
761 apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
762
763 map_molecule = mapping_atom_molecule(dcdr_env%lambda)
764 apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) &
765 = apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) + charge
766 END IF
767
768 ! And deallocate all the things!
769
770 CALL timestop(handle)
771 END SUBROUTINE apt_dr_localization
772
773END MODULE qs_dcdr
774
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:75
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
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
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:517
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
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:552
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1665
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:563
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:528
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:164
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:853
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:1063
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:731
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::molecule_of_atom
subroutine, public molecule_of_atom(molecule_set, atom_to_mol)
finds for each atom the molecule it belongs to
Definition molecule_types.F:423
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_dcdr_ao
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition qs_dcdr_ao.F:13
qs_dcdr_ao::core_dr
subroutine, public core_dr(qs_env, dcdr_env)
Core Hamiltonian contributions to the operator (the pseudopotentials)
Definition qs_dcdr_ao.F:399
qs_dcdr_ao::apply_op_constant_term
subroutine, public apply_op_constant_term(qs_env, dcdr_env, overlap1)
Build the perturbed density matrix correction depending on the overlap derivative.
Definition qs_dcdr_ao.F:117
qs_dcdr_ao::hr_mult_by_delta_1d
subroutine, public hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_or)
Enforce that one of the basis functions in < a | O | b > is centered on atom lambda.
Definition qs_dcdr_ao.F:703
qs_dcdr_ao::vhxc_r_perturbed_basis_functions
subroutine, public vhxc_r_perturbed_basis_functions(qs_env, dcdr_env)
The derivatives of the basis functions over which the HXC potential is integrated,...
Definition qs_dcdr_ao.F:632
qs_dcdr_ao::d_vhxc_dr
subroutine, public d_vhxc_dr(qs_env, dcdr_env)
The derivatives of the basis functions going into the HXC potential wrt nuclear positions.
Definition qs_dcdr_ao.F:504
qs_dcdr_ao::d_core_charge_density_dr
subroutine, public d_core_charge_density_dr(qs_env, dcdr_env)
Calculate the derivative of the Hartree term due to the core charge density.
Definition qs_dcdr_ao.F:327
qs_dcdr_utils
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition qs_dcdr_utils.F:13
qs_dcdr_utils::multiply_localization
subroutine, public multiply_localization(ao_matrix, mo_coeff, work, nmo, icenter, res)
Multiply (ao_matrix @ mo_coeff) and store the column icenter in res.
Definition qs_dcdr_utils.F:113
qs_dcdr_utils::dcdr_read_restart
subroutine, public dcdr_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_read_restart.
Definition qs_dcdr_utils.F:149
qs_dcdr_utils::shift_wannier_into_cell
subroutine, public shift_wannier_into_cell(r, cell, r_shifted)
...
Definition qs_dcdr_utils.F:527
qs_dcdr_utils::dcdr_write_restart
subroutine, public dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_write_restart.
Definition qs_dcdr_utils.F:299
qs_dcdr
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition qs_dcdr.F:13
qs_dcdr::prepare_per_atom
subroutine, public prepare_per_atom(dcdr_env, qs_env)
Prepare the environment for a choice of lambda.
Definition qs_dcdr.F:90
qs_dcdr::apt_dr_localization
subroutine, public apt_dr_localization(qs_env, dcdr_env)
Calculate atomic polar tensor using the localized dipole operator.
Definition qs_dcdr.F:519
qs_dcdr::apt_dr
subroutine, public apt_dr(qs_env, dcdr_env)
Calculate atomic polar tensor.
Definition qs_dcdr.F:376
qs_dcdr::dcdr_response_dr
subroutine, public dcdr_response_dr(dcdr_env, p_env, qs_env)
Get the dC/dR by solving the Sternheimer equation, using the op_dR matrix.
Definition qs_dcdr.F:269
qs_dcdr::dcdr_build_op_dr
subroutine, public dcdr_build_op_dr(dcdr_env, qs_env)
Build the operator for the position perturbation.
Definition qs_dcdr.F:194
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, zatom, 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_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_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
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_linres_types::get_polar_env
subroutine, public get_polar_env(polar_env, do_raman, do_periodic, dberry_psi0, polar, psi1_dberry, run_stopped)
...
Definition qs_linres_types.F:1242
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:560
qs_moments::dipole_deriv_ao
subroutine, public dipole_deriv_ao(qs_env, difdip, deltar, order, rcc)
...
Definition qs_moments.F:3094
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition qs_p_env_types.F:14
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
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:127
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_linres_types::dcdr_env_type
Definition qs_linres_types.F:247
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_linres_types::polar_env_type
Definition qs_linres_types.F:136
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
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