Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
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 : 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 :
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 72 : CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, 1, dcdr_env%ref_point)
180 :
181 72 : CALL timestop(handle)
182 72 : 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 216 : 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 216 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: opdr_sym
204 :
205 216 : CALL timeset(routineN, handle)
206 :
207 216 : nao = dcdr_env%nao
208 :
209 : ! allocate matrix for the sum of the perturbation terms of the operator (dbcsr matrix)
210 216 : NULLIFY (opdr_sym)
211 216 : CALL dbcsr_allocate_matrix_set(opdr_sym, 1)
212 216 : ALLOCATE (opdr_sym(1)%matrix)
213 216 : CALL dbcsr_copy(opdr_sym(1)%matrix, dcdr_env%matrix_s1(1)%matrix) ! symmetric
214 216 : CALL dbcsr_set(opdr_sym(1)%matrix, 0.0_dp)
215 :
216 468 : DO ispin = 1, dcdr_env%nspins
217 252 : nmo = dcdr_env%nmo(ispin)
218 :
219 252 : CALL apply_op_constant_term(qs_env, dcdr_env) ! dcdr_env%matrix_apply_op_constant
220 : ! Hartree and Exchange-Correlation contributions
221 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_core_charge_1(dcdr_env%beta)%matrix, zero, one)
222 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_d_vhxc_dR(dcdr_env%beta, ispin)%matrix, one, one)
223 252 : 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 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_hc(dcdr_env%beta)%matrix, one, one)
227 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_ppnl_1(dcdr_env%beta)%matrix, one, one)
228 252 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_apply_op_constant(ispin)%matrix, one, one)
229 :
230 252 : CALL dbcsr_desymmetrize(opdr_sym(1)%matrix, dcdr_env%hamiltonian1(1)%matrix)
231 252 : 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 252 : 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 252 : 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 252 : 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 252 : nmo, alpha=1.0_dp, beta=1.0_dp)
245 252 : CALL cp_fm_release(buf)
246 :
247 : ! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
248 252 : CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(ispin))
249 :
250 720 : IF (dcdr_env%z_matrix_method) THEN
251 54 : 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 216 : CALL dbcsr_deallocate_matrix_set(opdr_sym)
257 :
258 216 : CALL timestop(handle)
259 216 : 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 180 : 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 180 : 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 180 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
283 : TYPE(section_vals_type), POINTER :: lr_section
284 :
285 180 : CALL timeset(routineN, handle)
286 180 : NULLIFY (linres_control, lr_section, logger)
287 :
288 : CALL get_qs_env(qs_env=qs_env, &
289 : linres_control=linres_control, &
290 180 : mos=mos)
291 :
292 180 : logger => cp_get_default_logger()
293 180 : 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 180 : extension=".linresLog")
297 180 : IF (output_unit > 0) THEN
298 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
299 90 : "*** Self consistent optimization of the response wavefunction ***"
300 : END IF
301 :
302 : ! allocate the vectors
303 738 : ALLOCATE (psi0_order(dcdr_env%nspins))
304 738 : ALLOCATE (psi1(dcdr_env%nspins))
305 738 : ALLOCATE (h1_psi0(dcdr_env%nspins))
306 :
307 378 : DO ispin = 1, dcdr_env%nspins
308 198 : CALL cp_fm_create(psi1(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
309 198 : CALL cp_fm_create(h1_psi0(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
310 198 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
311 378 : psi0_order(ispin) = mo_coeff
312 : END DO
313 :
314 : ! Restart
315 180 : IF (linres_control%linres_restart) THEN
316 18 : CALL dcdr_read_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
317 : ELSE
318 342 : DO ispin = 1, dcdr_env%nspins
319 342 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
320 : END DO
321 : END IF
322 :
323 180 : IF (output_unit > 0) THEN
324 : WRITE (output_unit, "(T10,A,I4,A)") &
325 90 : "Response to the perturbation operator referring to atom ", dcdr_env%lambda, &
326 180 : " displaced in "//ACHAR(dcdr_env%beta + 119)
327 : END IF
328 378 : DO ispin = 1, dcdr_env%nspins
329 198 : CALL cp_fm_set_all(dcdr_env%dCR(ispin), 0.0_dp)
330 378 : CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), h1_psi0(ispin))
331 : END DO
332 :
333 180 : linres_control%lr_triplet = .FALSE. ! we do singlet response
334 180 : linres_control%do_kernel = .TRUE.
335 180 : 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 180 : output_unit, should_stop)
344 378 : DO ispin = 1, dcdr_env%nspins
345 378 : 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 180 : IF (linres_control%linres_restart) THEN
350 18 : CALL dcdr_write_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
351 : END IF
352 :
353 : ! clean up
354 378 : DO ispin = 1, dcdr_env%nspins
355 198 : CALL cp_fm_release(psi1(ispin))
356 378 : CALL cp_fm_release(h1_psi0(ispin))
357 : END DO
358 180 : DEALLOCATE (psi1, h1_psi0, psi0_order)
359 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
360 180 : "PRINT%PROGRAM_RUN_INFO")
361 :
362 180 : CALL timestop(handle)
363 :
364 540 : 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 540 : 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 180 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
387 : TYPE(cp_fm_type) :: overlap1_MO, tmp_fm_like_mos
388 180 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0, psi1_dBerry
389 : TYPE(cp_fm_type), POINTER :: mo_coeff
390 180 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
391 : TYPE(polar_env_type), POINTER :: polar_env
392 180 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
393 :
394 180 : apt_basis_derivative = 0._dp
395 180 : apt_coeff_derivative = 0._dp
396 :
397 180 : CALL timeset(routineN, handle)
398 :
399 180 : NULLIFY (qs_kind_set, particle_set)
400 : CALL get_qs_env(qs_env=qs_env, &
401 : qs_kind_set=qs_kind_set, &
402 180 : particle_set=particle_set)
403 :
404 180 : nao = dcdr_env%nao
405 180 : apt_el => dcdr_env%apt_el_dcdr
406 180 : apt_nuc => dcdr_env%apt_nuc_dcdr
407 :
408 180 : f_spin = 2._dp/dcdr_env%nspins
409 :
410 396 : DO ispin = 1, dcdr_env%nspins
411 : ! Compute S^(1,R)_(ij)
412 216 : CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
413 216 : CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
414 216 : nmo = dcdr_env%nmo(ispin)
415 216 : mo_coeff => dcdr_env%mo_coeff(ispin)
416 216 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
417 216 : 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 216 : 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 216 : 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 216 : -1.0_dp, dcdr_env%dCR_prime(ispin))
430 216 : CALL cp_fm_release(overlap1_MO)
431 :
432 864 : DO alpha = 1, 3
433 648 : IF (.NOT. dcdr_env%z_matrix_method) THEN
434 :
435 : ! FIRST CONTRIBUTION: dCR * moments * mo
436 486 : CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
437 486 : CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
438 486 : 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 486 : -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 486 : tmp_fm_like_mos, ncol=nmo)
444 486 : CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
445 :
446 486 : apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
447 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
448 486 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
449 : ELSE
450 162 : 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 162 : 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 162 : temp1)
460 :
461 : CALL cp_fm_trace(dcdr_env%matrix_m_alpha(dcdr_env%beta, ispin), &
462 : psi1_dBerry(alpha, ispin), &
463 162 : temp2)
464 :
465 162 : 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 162 : apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
472 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
473 162 : = 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 648 : 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 648 : tmp_fm_like_mos, ncol=nmo)
483 648 : 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 648 : apt_basis_derivative = -f_spin*apt_basis_derivative
487 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
488 864 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
489 :
490 : END DO ! alpha
491 :
492 612 : 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 180 : CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
497 180 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
498 180 : IF (.NOT. ghost) THEN
499 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
500 180 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
501 : END IF
502 :
503 : ! And deallocate all the things!
504 180 : CALL cp_fm_release(tmp_fm_like_mos)
505 180 : CALL cp_fm_release(overlap1_MO)
506 :
507 180 : CALL timestop(handle)
508 180 : 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 36 : 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 36 : 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 36 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements, diagonal_elements2
536 : REAL(dp), DIMENSION(3) :: distance, r_shifted
537 36 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
538 36 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
539 : TYPE(cell_type), POINTER :: cell
540 36 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
541 36 : 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 36 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
546 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
547 : TYPE(polar_env_type), POINTER :: polar_env
548 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
549 :
550 36 : CALL timeset(routineN, handle)
551 :
552 36 : 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 36 : cell=cell)
559 :
560 36 : nsubset = SIZE(molecule_set)
561 36 : natom = SIZE(particle_set)
562 36 : apt_el => dcdr_env%apt_el_dcdr
563 36 : apt_nuc => dcdr_env%apt_nuc_dcdr
564 36 : apt_subset => dcdr_env%apt_el_dcdr_per_subset
565 36 : apt_center => dcdr_env%apt_el_dcdr_per_center
566 :
567 : ! Map wannier functions to atoms
568 36 : IF (dcdr_env%nspins == 1) THEN
569 36 : max_nbr_center = dcdr_env%nbr_center(1)
570 : ELSE
571 0 : max_nbr_center = MAX(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
572 : END IF
573 144 : ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
574 108 : ALLOCATE (mapping_atom_molecule(natom))
575 36 : centers_set => dcdr_env%centers_set
576 :
577 72 : DO ispin = 1, dcdr_env%nspins
578 180 : 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 144 : r_shifted=r_shifted)
583 :
584 144 : smallest_r = HUGE(0._dp)
585 612 : DO i = 1, natom
586 432 : distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
587 1728 : tmp_r = SUM(distance**2)
588 576 : IF (tmp_r < smallest_r) THEN
589 216 : mapping_wannier_atom(icenter, ispin) = i
590 216 : smallest_r = tmp_r
591 : END IF
592 : END DO
593 : END DO
594 :
595 : ! Map atoms to molecules
596 36 : CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
597 72 : IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
598 20 : DO icenter = 1, dcdr_env%nbr_center(ispin)
599 16 : map_atom = mapping_wannier_atom(icenter, ispin)
600 20 : map_molecule = mapping_atom_molecule(map_atom)
601 : END DO
602 : END IF
603 : END DO !ispin
604 :
605 36 : nao = dcdr_env%nao
606 36 : f_spin = 2._dp/dcdr_env%nspins
607 :
608 72 : DO ispin = 1, dcdr_env%nspins
609 : ! Compute S^(1,R)_(ij)
610 :
611 36 : ALLOCATE (tmp_fm_like_mos)
612 36 : ALLOCATE (overlap1_MO)
613 36 : CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
614 36 : CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
615 36 : nmo = dcdr_env%nmo(ispin)
616 36 : mo_coeff => dcdr_env%mo_coeff(ispin)
617 36 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
618 36 : 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 36 : 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 36 : 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 36 : -1.0_dp, dcdr_env%dCR_prime(ispin))
631 36 : CALL cp_fm_release(overlap1_MO)
632 :
633 108 : ALLOCATE (diagonal_elements(nmo))
634 72 : ALLOCATE (diagonal_elements2(nmo))
635 :
636 : ! Allocate temporary matrices
637 36 : ALLOCATE (tmp_fm)
638 36 : ALLOCATE (tmp_fm_momo)
639 36 : ALLOCATE (tmp_fm_momo2)
640 36 : CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
641 36 : CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
642 36 : CALL cp_fm_create(tmp_fm_momo2, dcdr_env%momo_fm_struct(ispin)%struct)
643 :
644 : ! FIRST CONTRIBUTION: dCR * moments * mo
645 36 : this_factor = -2._dp*f_spin
646 144 : DO alpha = 1, 3
647 108 : IF (.NOT. dcdr_env%z_matrix_method) THEN
648 :
649 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
650 432 : CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
651 : CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
652 432 : 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 540 : 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 108 : 0.0_dp, tmp_fm_momo)
662 108 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
663 :
664 : ELSE
665 0 : 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 0 : 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.0_dp, tmp_fm_momo)
672 0 : 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 0 : psi1_dBerry(alpha, ispin), 0.0_dp, tmp_fm_momo2)
677 0 : CALL cp_fm_get_diag(tmp_fm_momo2, diagonal_elements2)
678 :
679 0 : diagonal_elements(:) = diagonal_elements(:) - diagonal_elements2(:)
680 : END IF
681 :
682 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
683 432 : map_atom = mapping_wannier_atom(icenter, ispin)
684 432 : map_molecule = mapping_atom_molecule(map_atom)
685 432 : tmp_aptcontr = this_factor*diagonal_elements(icenter)
686 :
687 : apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
688 432 : = 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 540 : = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
692 : END DO
693 :
694 540 : apt_coeff_derivative = this_factor*SUM(diagonal_elements)
695 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
696 144 : = 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 36 : this_factor = -f_spin
704 144 : DO alpha = 1, 3
705 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
706 : ! Build the AO matrix with the right wannier center as reference point
707 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(1, dcdr_env%beta)%matrix, 0._dp)
708 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(2, dcdr_env%beta)%matrix, 0._dp)
709 432 : 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 432 : 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 540 : 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 108 : 0.0_dp, tmp_fm_momo)
721 108 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
722 :
723 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
724 432 : map_atom = mapping_wannier_atom(icenter, ispin)
725 432 : map_molecule = mapping_atom_molecule(map_atom)
726 432 : tmp_aptcontr = this_factor*diagonal_elements(icenter)
727 :
728 : apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
729 432 : = 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 540 : = 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 540 : apt_basis_derivative = this_factor*SUM(diagonal_elements)
737 :
738 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
739 144 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
740 :
741 : END DO ! alpha
742 36 : DEALLOCATE (diagonal_elements)
743 36 : DEALLOCATE (diagonal_elements2)
744 :
745 36 : CALL cp_fm_release(tmp_fm)
746 36 : CALL cp_fm_release(tmp_fm_like_mos)
747 36 : CALL cp_fm_release(tmp_fm_momo)
748 36 : CALL cp_fm_release(tmp_fm_momo2)
749 36 : DEALLOCATE (overlap1_MO)
750 36 : DEALLOCATE (tmp_fm)
751 36 : DEALLOCATE (tmp_fm_like_mos)
752 36 : DEALLOCATE (tmp_fm_momo)
753 144 : DEALLOCATE (tmp_fm_momo2)
754 : END DO !ispin
755 :
756 : ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
757 36 : CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
758 36 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
759 36 : IF (.NOT. ghost) THEN
760 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
761 36 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
762 :
763 36 : map_molecule = mapping_atom_molecule(dcdr_env%lambda)
764 : apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) &
765 36 : = 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 36 : CALL timestop(handle)
771 108 : END SUBROUTINE apt_dR_localization
772 :
773 : END MODULE qs_dcdr
774 :
|
