Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate the derivatives of the MO coefficients wrt nuclear coordinates
10 : !> \author Sandra Luber, Edward Ditler
11 : ! **************************************************************************************************
12 :
13 : MODULE 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 : build_local_moments_der_matrix
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 :
81 : CONTAINS
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 72 : 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 72 : POINTER :: sab_all
98 72 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
99 72 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
100 :
101 72 : CALL timeset(routineN, handle)
102 :
103 72 : 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 72 : particle_set=particle_set)
108 :
109 72 : natom = SIZE(particle_set)
110 72 : IF (dcdr_env%distributed_origin) dcdr_env%ref_point(:) = particle_set(dcdr_env%lambda)%r(:)
111 :
112 936 : dcdr_env%delta_basis_function = 0._dp
113 288 : 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 288 : DO i = 1, 3
128 : ! S matrix
129 216 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
130 216 : CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_s1(1 + i)%matrix)
131 216 : 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 216 : 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 216 : sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
137 :
138 216 : CALL dbcsr_add(dcdr_env%matrix_s1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
139 216 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
140 :
141 : ! T matrix
142 216 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
143 216 : CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_t1(1 + i)%matrix)
144 216 : 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 216 : 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 216 : sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
150 :
151 216 : CALL dbcsr_add(dcdr_env%matrix_t1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
152 288 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
153 : END DO
154 :
155 : ! Operator:
156 156 : DO ispin = 1, dcdr_env%nspins
157 408 : DO i = 1, 3
158 252 : CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
159 252 : CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
160 252 : CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i)%matrix, 0.0_dp)
161 252 : CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i + 3)%matrix, 0.0_dp)
162 252 : CALL dbcsr_set(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, 0.0_dp)
163 336 : CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
164 : END DO
165 : END DO
166 :
167 72 : CALL core_dR(qs_env, dcdr_env) ! dcdr_env%matrix_ppnl_1, hc
168 72 : CALL d_vhxc_dR(qs_env, dcdr_env) ! dcdr_env%matrix_d_vhxc_dR
169 72 : CALL d_core_charge_density_dR(qs_env, dcdr_env) ! dcdr_env%matrix_core_charge_1
170 72 : CALL vhxc_R_perturbed_basis_functions(qs_env, dcdr_env) ! dcdr_env%matrix_vhxc_perturbed_basis
171 :
172 : ! APT:
173 288 : DO i = 1, 3
174 936 : DO j = 1, 3
175 864 : CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0._dp)
176 : END DO
177 : END DO
178 :
179 : CALL build_local_moments_der_matrix(qs_env, dcdr_env%matrix_difdip, 1, 0, &
180 : ref_point=dcdr_env%ref_point, ordered=.TRUE., &
181 72 : deltaR=dcdr_env%delta_basis_function, neighbor_image=.FALSE.)
182 :
183 72 : CALL timestop(handle)
184 72 : END SUBROUTINE prepare_per_atom
185 :
186 : ! **************************************************************************************************
187 : !> \brief Build the operator for the position perturbation
188 : !> \param dcdr_env ...
189 : !> \param qs_env ...
190 : !> \authors Sandra Luber
191 : !> Edward Ditler
192 : !> Ravi Kumar
193 : !> Rangsiman Ketkaew
194 : ! **************************************************************************************************
195 216 : SUBROUTINE dcdr_build_op_dR(dcdr_env, qs_env)
196 :
197 : TYPE(dcdr_env_type) :: dcdr_env
198 : TYPE(qs_environment_type), POINTER :: qs_env
199 :
200 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_build_op_dR'
201 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
202 :
203 : INTEGER :: handle, ispin, nao, nmo
204 : TYPE(cp_fm_type) :: buf
205 216 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: opdr_sym
206 :
207 216 : CALL timeset(routineN, handle)
208 :
209 216 : nao = dcdr_env%nao
210 :
211 : ! allocate matrix for the sum of the perturbation terms of the operator (dbcsr matrix)
212 216 : NULLIFY (opdr_sym)
213 216 : CALL dbcsr_allocate_matrix_set(opdr_sym, 1)
214 216 : ALLOCATE (opdr_sym(1)%matrix)
215 216 : CALL dbcsr_copy(opdr_sym(1)%matrix, dcdr_env%matrix_s1(1)%matrix) ! symmetric
216 216 : CALL dbcsr_set(opdr_sym(1)%matrix, 0.0_dp)
217 :
218 468 : DO ispin = 1, dcdr_env%nspins
219 252 : nmo = dcdr_env%nmo(ispin)
220 :
221 252 : CALL apply_op_constant_term(qs_env, dcdr_env) ! dcdr_env%matrix_apply_op_constant
222 : ! Hartree and Exchange-Correlation contributions
223 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_core_charge_1(dcdr_env%beta)%matrix, zero, one)
224 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_d_vhxc_dR(dcdr_env%beta, ispin)%matrix, one, one)
225 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_vhxc_perturbed_basis(ispin, dcdr_env%beta)%matrix, one, one)
226 :
227 : ! Core Hamiltonian contributions
228 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_hc(dcdr_env%beta)%matrix, one, one)
229 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_ppnl_1(dcdr_env%beta)%matrix, one, one)
230 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_apply_op_constant(ispin)%matrix, one, one)
231 :
232 252 : CALL dbcsr_desymmetrize(opdr_sym(1)%matrix, dcdr_env%hamiltonian1(1)%matrix)
233 252 : CALL dbcsr_add(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%matrix_t1(dcdr_env%beta + 1)%matrix, one, one)
234 :
235 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%mo_coeff(ispin), &
236 252 : dcdr_env%op_dR(ispin), ncol=nmo)
237 :
238 : ! The overlap derivative terms for the Sternheimer equation
239 : ! buf = mo * (-mo * matrix_ks * mo)
240 252 : CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
241 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
242 : -1.0_dp, dcdr_env%mo_coeff(ispin), dcdr_env%chc(ispin), &
243 252 : 0.0_dp, buf)
244 :
245 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, buf, dcdr_env%op_dR(ispin), &
246 252 : nmo, alpha=1.0_dp, beta=1.0_dp)
247 252 : CALL cp_fm_release(buf)
248 :
249 : ! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
250 252 : CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(ispin))
251 :
252 720 : IF (dcdr_env%z_matrix_method) THEN
253 54 : CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin))
254 : END IF
255 :
256 : END DO
257 :
258 216 : CALL dbcsr_deallocate_matrix_set(opdr_sym)
259 :
260 216 : CALL timestop(handle)
261 216 : END SUBROUTINE dcdr_build_op_dR
262 :
263 : ! **************************************************************************************************
264 : !> \brief Get the dC/dR by solving the Sternheimer equation, using the op_dR matrix
265 : !> \param dcdr_env ...
266 : !> \param p_env ...
267 : !> \param qs_env ...
268 : !> \authors SL, ED
269 : ! **************************************************************************************************
270 180 : SUBROUTINE dcdr_response_dR(dcdr_env, p_env, qs_env)
271 :
272 : TYPE(dcdr_env_type) :: dcdr_env
273 : TYPE(qs_p_env_type) :: p_env
274 : TYPE(qs_environment_type), POINTER :: qs_env
275 :
276 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_response_dR'
277 :
278 : INTEGER :: handle, ispin, output_unit
279 : LOGICAL :: should_stop
280 180 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
281 : TYPE(cp_fm_type), POINTER :: mo_coeff
282 : TYPE(cp_logger_type), POINTER :: logger
283 : TYPE(linres_control_type), POINTER :: linres_control
284 180 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
285 : TYPE(section_vals_type), POINTER :: lr_section
286 :
287 180 : CALL timeset(routineN, handle)
288 180 : NULLIFY (linres_control, lr_section, logger)
289 :
290 : CALL get_qs_env(qs_env=qs_env, &
291 : linres_control=linres_control, &
292 180 : mos=mos)
293 :
294 180 : logger => cp_get_default_logger()
295 180 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
296 :
297 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
298 180 : extension=".linresLog")
299 180 : IF (output_unit > 0) THEN
300 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
301 90 : "*** Self consistent optimization of the response wavefunction ***"
302 : END IF
303 :
304 : ! allocate the vectors
305 738 : ALLOCATE (psi0_order(dcdr_env%nspins))
306 558 : ALLOCATE (psi1(dcdr_env%nspins))
307 558 : ALLOCATE (h1_psi0(dcdr_env%nspins))
308 :
309 378 : DO ispin = 1, dcdr_env%nspins
310 198 : CALL cp_fm_create(psi1(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
311 198 : CALL cp_fm_create(h1_psi0(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
312 198 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
313 378 : psi0_order(ispin) = mo_coeff
314 : END DO
315 :
316 378 : DO ispin = 1, dcdr_env%nspins
317 198 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
318 378 : CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
319 : END DO
320 : ! Restart
321 180 : IF (linres_control%linres_restart) THEN
322 18 : CALL dcdr_read_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
323 : END IF
324 :
325 180 : IF (output_unit > 0) THEN
326 : WRITE (output_unit, "(T10,A,I4,A)") &
327 90 : "Response to the perturbation operator referring to atom ", dcdr_env%lambda, &
328 180 : " displaced in "//ACHAR(dcdr_env%beta + 119)
329 : END IF
330 378 : DO ispin = 1, dcdr_env%nspins
331 198 : CALL cp_fm_set_all(dcdr_env%dCR(ispin), 0.0_dp)
332 378 : CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), h1_psi0(ispin))
333 : END DO
334 :
335 180 : linres_control%lr_triplet = .FALSE. ! we do singlet response
336 180 : linres_control%do_kernel = .TRUE.
337 180 : linres_control%converged = .FALSE.
338 :
339 : ! Position perturbation to get dCR
340 : ! (H0-E0) psi1 = (H1-E1) psi0
341 : ! psi1 = the perturbed wavefunction
342 : ! h1_psi0 = (H1-E1-S1*\varepsilon)
343 : ! psi0_order = the unperturbed wavefunction
344 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
345 180 : output_unit, should_stop)
346 378 : DO ispin = 1, dcdr_env%nspins
347 378 : CALL cp_fm_to_fm(psi1(ispin), dcdr_env%dCR(ispin))
348 : END DO
349 :
350 : ! Write the new result to the restart file
351 180 : IF (linres_control%linres_restart) THEN
352 18 : CALL dcdr_write_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
353 : END IF
354 :
355 : ! clean up
356 378 : DO ispin = 1, dcdr_env%nspins
357 198 : CALL cp_fm_release(psi1(ispin))
358 378 : CALL cp_fm_release(h1_psi0(ispin))
359 : END DO
360 180 : DEALLOCATE (psi1, h1_psi0, psi0_order)
361 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
362 180 : "PRINT%PROGRAM_RUN_INFO")
363 :
364 180 : CALL timestop(handle)
365 :
366 540 : END SUBROUTINE dcdr_response_dR
367 :
368 : ! **************************************************************************************************
369 : !> \brief Calculate atomic polar tensor
370 : !> \param qs_env ...
371 : !> \param dcdr_env ...
372 : !> \authors Sandra Luber
373 : !> Edward Ditler
374 : !> Ravi Kumar
375 : !> Rangsiman Ketkaew
376 : ! **************************************************************************************************
377 540 : SUBROUTINE apt_dR(qs_env, dcdr_env)
378 : TYPE(qs_environment_type), POINTER :: qs_env
379 : TYPE(dcdr_env_type) :: dcdr_env
380 :
381 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR'
382 :
383 : INTEGER :: alpha, handle, ikind, ispin, nao, nmo
384 : LOGICAL :: ghost
385 : REAL(dp) :: apt_basis_derivative, &
386 : apt_coeff_derivative, charge, f_spin, &
387 : temp1, temp2
388 180 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
389 : TYPE(cp_fm_type) :: overlap1_MO, tmp_fm_like_mos
390 180 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
391 : TYPE(cp_fm_type), POINTER :: mo_coeff
392 180 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
393 : TYPE(polar_env_type), POINTER :: polar_env
394 180 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
395 :
396 : apt_basis_derivative = 0._dp
397 : apt_coeff_derivative = 0._dp
398 :
399 180 : CALL timeset(routineN, handle)
400 :
401 180 : NULLIFY (qs_kind_set, particle_set)
402 : CALL get_qs_env(qs_env=qs_env, &
403 : qs_kind_set=qs_kind_set, &
404 180 : particle_set=particle_set)
405 :
406 180 : nao = dcdr_env%nao
407 180 : apt_el => dcdr_env%apt_el_dcdr
408 180 : apt_nuc => dcdr_env%apt_nuc_dcdr
409 :
410 180 : f_spin = 2._dp/dcdr_env%nspins
411 :
412 396 : DO ispin = 1, dcdr_env%nspins
413 : ! Compute S^(1,R)_(ij)
414 216 : CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
415 216 : CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
416 216 : nmo = dcdr_env%nmo(ispin)
417 216 : mo_coeff => dcdr_env%mo_coeff(ispin)
418 216 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
419 216 : CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
420 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
421 216 : tmp_fm_like_mos, ncol=nmo)
422 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
423 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
424 216 : 0.0_dp, overlap1_MO)
425 :
426 : ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
427 : ! We get the negative of the coefficients out of the linres solver
428 : ! And apply the constant correction due to the overlap derivative.
429 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
430 : -0.5_dp, mo_coeff, overlap1_MO, &
431 216 : -1.0_dp, dcdr_env%dCR_prime(ispin))
432 216 : CALL cp_fm_release(overlap1_MO)
433 :
434 864 : DO alpha = 1, 3
435 648 : IF (.NOT. dcdr_env%z_matrix_method) THEN
436 :
437 : ! FIRST CONTRIBUTION: dCR * moments * mo
438 486 : CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
439 486 : CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
440 486 : CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
441 : CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
442 486 : -dcdr_env%ref_point(alpha), 1._dp)
443 :
444 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(ispin), &
445 486 : tmp_fm_like_mos, ncol=nmo)
446 486 : CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
447 :
448 486 : apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
449 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
450 486 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
451 : ELSE
452 162 : CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
453 : CALL get_polar_env(polar_env=polar_env, psi1_dBerry=psi1_dBerry, &
454 162 : dBerry_psi0=dBerry_psi0)
455 :
456 : ! Note that here dcdr_env%dCR_prime contains only occ-occ block contribution,
457 : ! dcdr_env%dCR(ispin) is zero because we didn't run response calculation for dcdR.
458 :
459 : CALL cp_fm_trace(dBerry_psi0(alpha, ispin), &
460 : dcdr_env%dCR_prime(ispin), &
461 162 : temp1)
462 :
463 : CALL cp_fm_trace(dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
464 : psi1_dBerry(alpha, ispin), &
465 162 : temp2)
466 :
467 162 : apt_coeff_derivative = temp1 - temp2
468 :
469 : ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
470 : ! - apt_coeff_derivative , here the trace is negative to compensate the
471 : ! -ve sign in APTs= - 2 Z. M_alpha
472 :
473 162 : apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
474 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
475 162 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
476 : END IF
477 :
478 : ! SECOND CONTRIBUTION: We assemble all combinations of r_i, d(chi)/d(idir)
479 : ! difdip contains derivatives with respect to atom dcdr_env%lambda
480 : ! difdip(alpha, beta): < a | r_alpha | db/dR_beta >
481 : ! Multiply by the MO coefficients
482 648 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
483 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, mo_coeff, &
484 648 : tmp_fm_like_mos, ncol=nmo)
485 648 : CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_basis_derivative)
486 :
487 : ! The negative sign compensates for the nuclear-coordinate derivative convention.
488 648 : apt_basis_derivative = -f_spin*apt_basis_derivative
489 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
490 864 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
491 :
492 : END DO ! alpha
493 :
494 612 : CALL cp_fm_release(tmp_fm_like_mos)
495 : END DO !ispin
496 :
497 : ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
498 180 : CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
499 180 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
500 180 : IF (.NOT. ghost) THEN
501 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
502 180 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
503 : END IF
504 :
505 : ! And deallocate all the things!
506 180 : CALL cp_fm_release(tmp_fm_like_mos)
507 180 : CALL cp_fm_release(overlap1_MO)
508 :
509 180 : CALL timestop(handle)
510 180 : END SUBROUTINE apt_dR
511 :
512 : ! **************************************************************************************************
513 : !> \brief Calculate atomic polar tensor using the localized dipole operator
514 : !> \param qs_env ...
515 : !> \param dcdr_env ...
516 : !> \authors Edward Ditler
517 : !> Ravi Kumar
518 : !> Rangsiman Ketkaew
519 : ! **************************************************************************************************
520 36 : SUBROUTINE apt_dR_localization(qs_env, dcdr_env)
521 : TYPE(qs_environment_type), POINTER :: qs_env
522 : TYPE(dcdr_env_type) :: dcdr_env
523 :
524 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR_localization'
525 :
526 : INTEGER :: alpha, handle, i, icenter, ikind, ispin, &
527 : map_atom, map_molecule, &
528 : max_nbr_center, nao, natom, nmo, &
529 : nsubset
530 : INTEGER, ALLOCATABLE, DIMENSION(:) :: mapping_atom_molecule
531 36 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mapping_wannier_atom
532 : LOGICAL :: ghost
533 : REAL(dp) :: apt_basis_derivative, &
534 : apt_coeff_derivative, charge, f_spin, &
535 : smallest_r, this_factor, tmp_aptcontr, &
536 : tmp_r
537 36 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements, diagonal_elements2
538 : REAL(dp), DIMENSION(3) :: distance, r_shifted
539 36 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
540 36 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
541 : TYPE(cell_type), POINTER :: cell
542 36 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
543 36 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
544 : TYPE(cp_fm_type), POINTER :: mo_coeff, overlap1_MO, tmp_fm, &
545 : tmp_fm_like_mos, tmp_fm_momo, &
546 : tmp_fm_momo2
547 36 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
548 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
549 : TYPE(polar_env_type), POINTER :: polar_env
550 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
551 :
552 36 : CALL timeset(routineN, handle)
553 :
554 36 : NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
555 :
556 : CALL get_qs_env(qs_env=qs_env, &
557 : qs_kind_set=qs_kind_set, &
558 : particle_set=particle_set, &
559 : molecule_set=molecule_set, &
560 36 : cell=cell)
561 :
562 36 : nsubset = SIZE(molecule_set)
563 36 : natom = SIZE(particle_set)
564 36 : apt_el => dcdr_env%apt_el_dcdr
565 36 : apt_nuc => dcdr_env%apt_nuc_dcdr
566 36 : apt_subset => dcdr_env%apt_el_dcdr_per_subset
567 36 : apt_center => dcdr_env%apt_el_dcdr_per_center
568 :
569 : ! Map wannier functions to atoms
570 36 : IF (dcdr_env%nspins == 1) THEN
571 36 : max_nbr_center = dcdr_env%nbr_center(1)
572 : ELSE
573 0 : max_nbr_center = MAX(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
574 : END IF
575 144 : ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
576 108 : ALLOCATE (mapping_atom_molecule(natom))
577 36 : centers_set => dcdr_env%centers_set
578 :
579 72 : DO ispin = 1, dcdr_env%nspins
580 180 : DO icenter = 1, dcdr_env%nbr_center(ispin)
581 : ! For every center we check which atom is closest
582 : CALL shift_wannier_into_cell(r=centers_set(ispin)%array(1:3, icenter), &
583 : cell=cell, &
584 144 : r_shifted=r_shifted)
585 :
586 144 : smallest_r = HUGE(0._dp)
587 612 : DO i = 1, natom
588 432 : distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
589 1728 : tmp_r = SUM(distance**2)
590 576 : IF (tmp_r < smallest_r) THEN
591 216 : mapping_wannier_atom(icenter, ispin) = i
592 216 : smallest_r = tmp_r
593 : END IF
594 : END DO
595 : END DO
596 :
597 : ! Map atoms to molecules
598 36 : CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
599 72 : IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
600 20 : DO icenter = 1, dcdr_env%nbr_center(ispin)
601 16 : map_atom = mapping_wannier_atom(icenter, ispin)
602 20 : map_molecule = mapping_atom_molecule(map_atom)
603 : END DO
604 : END IF
605 : END DO !ispin
606 :
607 36 : nao = dcdr_env%nao
608 36 : f_spin = 2._dp/dcdr_env%nspins
609 :
610 72 : DO ispin = 1, dcdr_env%nspins
611 : ! Compute S^(1,R)_(ij)
612 :
613 36 : ALLOCATE (tmp_fm_like_mos)
614 36 : ALLOCATE (overlap1_MO)
615 36 : CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
616 36 : CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
617 36 : nmo = dcdr_env%nmo(ispin)
618 36 : mo_coeff => dcdr_env%mo_coeff(ispin)
619 36 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
620 36 : CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
621 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
622 36 : tmp_fm_like_mos, ncol=nmo)
623 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
624 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
625 36 : 0.0_dp, overlap1_MO)
626 :
627 : ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
628 : ! We get the negative of the coefficients out of the linres solver
629 : ! And apply the constant correction due to the overlap derivative.
630 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
631 : -0.5_dp, mo_coeff, overlap1_MO, &
632 36 : -1.0_dp, dcdr_env%dCR_prime(ispin))
633 36 : CALL cp_fm_release(overlap1_MO)
634 :
635 108 : ALLOCATE (diagonal_elements(nmo))
636 72 : ALLOCATE (diagonal_elements2(nmo))
637 :
638 : ! Allocate temporary matrices
639 36 : ALLOCATE (tmp_fm)
640 36 : ALLOCATE (tmp_fm_momo)
641 36 : ALLOCATE (tmp_fm_momo2)
642 36 : CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
643 36 : CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
644 36 : CALL cp_fm_create(tmp_fm_momo2, dcdr_env%momo_fm_struct(ispin)%struct)
645 :
646 : ! FIRST CONTRIBUTION: dCR * moments * mo
647 36 : this_factor = -2._dp*f_spin
648 144 : DO alpha = 1, 3
649 108 : IF (.NOT. dcdr_env%z_matrix_method) THEN
650 :
651 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
652 432 : CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
653 : CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
654 432 : ref_point=centers_set(ispin)%array(1:3, icenter))
655 : CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
656 : mo_coeff=dcdr_env%dCR_prime(ispin), work=tmp_fm, nmo=nmo, &
657 : icenter=icenter, &
658 540 : res=tmp_fm_like_mos)
659 : END DO
660 :
661 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
662 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
663 108 : 0.0_dp, tmp_fm_momo)
664 108 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
665 :
666 : ELSE
667 0 : CALL get_qs_env(qs_env=qs_env, polar_env=polar_env)
668 : CALL get_polar_env(polar_env=polar_env, psi1_dBerry=psi1_dBerry, &
669 0 : dBerry_psi0=dBerry_psi0)
670 :
671 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
672 : 1.0_dp, dcdr_env%dCR_prime(ispin), dBerry_psi0(alpha, ispin), &
673 0 : 0.0_dp, tmp_fm_momo)
674 0 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
675 :
676 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
677 : 1.0_dp, dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
678 0 : psi1_dBerry(alpha, ispin), 0.0_dp, tmp_fm_momo2)
679 0 : CALL cp_fm_get_diag(tmp_fm_momo2, diagonal_elements2)
680 :
681 0 : diagonal_elements(:) = diagonal_elements(:) - diagonal_elements2(:)
682 : END IF
683 :
684 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
685 432 : map_atom = mapping_wannier_atom(icenter, ispin)
686 432 : map_molecule = mapping_atom_molecule(map_atom)
687 432 : tmp_aptcontr = this_factor*diagonal_elements(icenter)
688 :
689 : apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
690 432 : = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
691 :
692 : apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
693 540 : = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
694 : END DO
695 :
696 540 : apt_coeff_derivative = this_factor*SUM(diagonal_elements)
697 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
698 144 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
699 : END DO
700 :
701 : ! SECOND CONTRIBUTION: We assemble all combinations of r_i, dphi/d(idir)
702 : ! build part with AOs differentiated with respect to nuclear coordinates
703 : ! difdip contains derivatives with respect to atom dcdr_env%lambda
704 : ! difdip(alpha, beta): < a | r_alpha | d b/dR_beta >
705 36 : this_factor = -f_spin
706 144 : DO alpha = 1, 3
707 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
708 : ! Build the AO matrix with the right wannier center as reference point
709 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(1, dcdr_env%beta)%matrix, 0._dp)
710 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(2, dcdr_env%beta)%matrix, 0._dp)
711 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(3, dcdr_env%beta)%matrix, 0._dp)
712 : CALL build_local_moments_der_matrix(qs_env, dcdr_env%matrix_difdip, 1, 0, &
713 : ref_point=centers_set(ispin)%array(1:3, icenter), &
714 : ordered=.TRUE., deltaR=dcdr_env%delta_basis_function, &
715 432 : neighbor_image=.FALSE.)
716 : CALL multiply_localization(ao_matrix=dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, &
717 : mo_coeff=mo_coeff, work=tmp_fm, nmo=nmo, &
718 : icenter=icenter, &
719 540 : res=tmp_fm_like_mos)
720 : END DO ! icenter
721 :
722 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
723 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
724 108 : 0.0_dp, tmp_fm_momo)
725 108 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
726 :
727 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
728 432 : map_atom = mapping_wannier_atom(icenter, ispin)
729 432 : map_molecule = mapping_atom_molecule(map_atom)
730 432 : tmp_aptcontr = this_factor*diagonal_elements(icenter)
731 :
732 : apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
733 432 : = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
734 :
735 : apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
736 540 : = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
737 : END DO
738 :
739 : ! The negative sign compensates for the nuclear-coordinate derivative convention.
740 540 : apt_basis_derivative = this_factor*SUM(diagonal_elements)
741 :
742 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
743 144 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
744 :
745 : END DO ! alpha
746 36 : DEALLOCATE (diagonal_elements)
747 36 : DEALLOCATE (diagonal_elements2)
748 :
749 36 : CALL cp_fm_release(tmp_fm)
750 36 : CALL cp_fm_release(tmp_fm_like_mos)
751 36 : CALL cp_fm_release(tmp_fm_momo)
752 36 : CALL cp_fm_release(tmp_fm_momo2)
753 36 : DEALLOCATE (overlap1_MO)
754 36 : DEALLOCATE (tmp_fm)
755 36 : DEALLOCATE (tmp_fm_like_mos)
756 36 : DEALLOCATE (tmp_fm_momo)
757 144 : DEALLOCATE (tmp_fm_momo2)
758 : END DO !ispin
759 :
760 : ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
761 36 : CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
762 36 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
763 36 : IF (.NOT. ghost) THEN
764 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
765 36 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
766 :
767 36 : map_molecule = mapping_atom_molecule(dcdr_env%lambda)
768 : apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) &
769 36 : = apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) + charge
770 : END IF
771 :
772 : ! And deallocate all the things!
773 :
774 36 : CALL timestop(handle)
775 108 : END SUBROUTINE apt_dR_localization
776 :
777 : END MODULE qs_dcdr
|
