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_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
20 : dbcsr_allocate_matrix_set,&
21 : dbcsr_deallocate_matrix_set
22 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
23 : cp_fm_scale_and_add,&
24 : cp_fm_trace
25 : USE cp_fm_types, ONLY: cp_fm_create,&
26 : cp_fm_get_diag,&
27 : cp_fm_release,&
28 : cp_fm_set_all,&
29 : cp_fm_to_fm,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
34 : cp_print_key_unit_nr
35 : USE dbcsr_api, ONLY: dbcsr_add,&
36 : dbcsr_copy,&
37 : dbcsr_desymmetrize,&
38 : dbcsr_p_type,&
39 : dbcsr_set
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 : linres_control_type
64 : USE qs_mo_types, ONLY: get_mo_set,&
65 : mo_set_type
66 : USE qs_moments, ONLY: build_local_moment_matrix,&
67 : dipole_deriv_ao
68 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
69 : USE qs_p_env_types, ONLY: qs_p_env_type
70 : #include "./base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 :
74 : PRIVATE
75 : PUBLIC :: prepare_per_atom, dcdr_response_dR, dcdr_build_op_dR, apt_dR, apt_dR_localization
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr'
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief Prepare the environment for a choice of lambda
83 : !> \param dcdr_env ...
84 : !> \param qs_env ...
85 : !> \author Edward Ditler
86 : ! **************************************************************************************************
87 36 : SUBROUTINE prepare_per_atom(dcdr_env, qs_env)
88 : TYPE(dcdr_env_type) :: dcdr_env
89 : TYPE(qs_environment_type), POINTER :: qs_env
90 :
91 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_per_atom'
92 :
93 : INTEGER :: handle, i, ispin, j, natom
94 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
95 36 : POINTER :: sab_all
96 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
97 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
98 :
99 36 : CALL timeset(routineN, handle)
100 :
101 36 : NULLIFY (sab_all, qs_kind_set, particle_set)
102 : CALL get_qs_env(qs_env=qs_env, &
103 : sab_all=sab_all, &
104 : qs_kind_set=qs_kind_set, &
105 36 : particle_set=particle_set)
106 :
107 36 : natom = SIZE(particle_set)
108 36 : IF (dcdr_env%distributed_origin) dcdr_env%ref_point(:) = particle_set(dcdr_env%lambda)%r(:)
109 :
110 468 : dcdr_env%delta_basis_function = 0._dp
111 144 : dcdr_env%delta_basis_function(:, dcdr_env%lambda) = 1._dp
112 :
113 : ! S matrix
114 : ! S1 = - < da/dr | b > * delta_a - < a | db/dr > * delta_b
115 :
116 : ! matrix_s(2:4) are anti-symmetric matrices and contain derivatives wrt. to < a |
117 : ! = < da/dR | b > = - < da/dr | b > = < a | db/dr >
118 : ! matrix_s1(2:4) = d/dR < a | b >
119 : ! and it's built as
120 : ! = - matrix_s * delta_b + matrix_s * delta_a
121 : ! = - < da/dR | b > * delta_b + < da/dR | b > * delta_a
122 : ! = + < da/dr | b > * delta_b - < da/dr | b > * delta_a
123 : ! = - < a | db/dr > * delta_b - < da/dr | b > * delta_a
124 :
125 144 : DO i = 1, 3
126 : ! S matrix
127 108 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
128 108 : CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_s1(1 + i)%matrix)
129 108 : CALL dbcsr_desymmetrize(dcdr_env%matrix_s(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
130 :
131 : CALL hr_mult_by_delta_1d(dcdr_env%matrix_s1(1 + i)%matrix, qs_kind_set, "ORB", &
132 108 : sab_all, dcdr_env%lambda, direction_Or=.TRUE.)
133 : CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
134 108 : sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
135 :
136 108 : CALL dbcsr_add(dcdr_env%matrix_s1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
137 108 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
138 :
139 : ! T matrix
140 108 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
141 108 : CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_t1(1 + i)%matrix)
142 108 : CALL dbcsr_desymmetrize(dcdr_env%matrix_t(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix)
143 :
144 : CALL hr_mult_by_delta_1d(dcdr_env%matrix_t1(1 + i)%matrix, qs_kind_set, "ORB", &
145 108 : sab_all, dcdr_env%lambda, direction_Or=.TRUE.)
146 : CALL hr_mult_by_delta_1d(dcdr_env%matrix_nosym_temp(i)%matrix, qs_kind_set, "ORB", &
147 108 : sab_all, dcdr_env%lambda, direction_Or=.FALSE.)
148 :
149 108 : CALL dbcsr_add(dcdr_env%matrix_t1(1 + i)%matrix, dcdr_env%matrix_nosym_temp(i)%matrix, -1._dp, +1._dp)
150 144 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
151 : END DO
152 :
153 : ! Operator:
154 78 : DO ispin = 1, dcdr_env%nspins
155 204 : DO i = 1, 3
156 126 : CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
157 126 : CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
158 126 : CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i)%matrix, 0.0_dp)
159 126 : CALL dbcsr_set(dcdr_env%matrix_vhxc_perturbed_basis(ispin, i + 3)%matrix, 0.0_dp)
160 126 : CALL dbcsr_set(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, 0.0_dp)
161 168 : CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
162 : END DO
163 : END DO
164 :
165 36 : CALL core_dR(qs_env, dcdr_env) ! dcdr_env%matrix_ppnl_1, hc
166 36 : CALL d_vhxc_dR(qs_env, dcdr_env) ! dcdr_env%matrix_d_vhxc_dR
167 36 : CALL d_core_charge_density_dR(qs_env, dcdr_env) ! dcdr_env%matrix_core_charge_1
168 36 : CALL vhxc_R_perturbed_basis_functions(qs_env, dcdr_env) ! dcdr_env%matrix_vhxc_perturbed_basis
169 :
170 : ! APT:
171 144 : DO i = 1, 3
172 468 : DO j = 1, 3
173 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0._dp)
174 : END DO
175 : END DO
176 :
177 36 : CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, 1, dcdr_env%ref_point)
178 :
179 36 : CALL timestop(handle)
180 36 : END SUBROUTINE prepare_per_atom
181 :
182 : ! **************************************************************************************************
183 : !> \brief Build the operator for the position perturbation
184 : !> \param dcdr_env ...
185 : !> \param qs_env ...
186 : !> \authors SL, ED
187 : ! **************************************************************************************************
188 108 : SUBROUTINE dcdr_build_op_dR(dcdr_env, qs_env)
189 :
190 : TYPE(dcdr_env_type) :: dcdr_env
191 : TYPE(qs_environment_type), POINTER :: qs_env
192 :
193 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_build_op_dR'
194 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
195 :
196 : INTEGER :: handle, ispin, nao, nmo
197 : TYPE(cp_fm_type) :: buf
198 108 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: opdr_sym
199 :
200 108 : CALL timeset(routineN, handle)
201 :
202 108 : nao = dcdr_env%nao
203 :
204 : ! allocate matrix for the sum of the perturbation terms of the operator (dbcsr matrix)
205 108 : NULLIFY (opdr_sym)
206 108 : CALL dbcsr_allocate_matrix_set(opdr_sym, 1)
207 108 : ALLOCATE (opdr_sym(1)%matrix)
208 108 : CALL dbcsr_copy(opdr_sym(1)%matrix, dcdr_env%matrix_s1(1)%matrix) ! symmetric
209 108 : CALL dbcsr_set(opdr_sym(1)%matrix, 0.0_dp)
210 :
211 234 : DO ispin = 1, dcdr_env%nspins
212 126 : nmo = dcdr_env%nmo(ispin)
213 :
214 126 : CALL apply_op_constant_term(qs_env, dcdr_env) ! dcdr_env%matrix_apply_op_constant
215 : ! Hartree and Exchange-Correlation contributions
216 126 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_core_charge_1(dcdr_env%beta)%matrix, zero, one)
217 126 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_d_vhxc_dR(dcdr_env%beta, ispin)%matrix, one, one)
218 126 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_vhxc_perturbed_basis(ispin, dcdr_env%beta)%matrix, one, one)
219 :
220 : ! Core Hamiltonian contributions
221 126 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_hc(dcdr_env%beta)%matrix, one, one)
222 126 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_ppnl_1(dcdr_env%beta)%matrix, one, one)
223 126 : CALL dbcsr_add(opdr_sym(1)%matrix, dcdr_env%matrix_apply_op_constant(ispin)%matrix, one, one)
224 :
225 126 : CALL dbcsr_desymmetrize(opdr_sym(1)%matrix, dcdr_env%hamiltonian1(1)%matrix)
226 126 : CALL dbcsr_add(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%matrix_t1(dcdr_env%beta + 1)%matrix, one, one)
227 :
228 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%hamiltonian1(1)%matrix, dcdr_env%mo_coeff(ispin), &
229 126 : dcdr_env%op_dR(ispin), ncol=nmo)
230 :
231 : ! The overlap derivative terms for the Sternheimer equation
232 : ! buf = mo * (-mo * matrix_ks * mo)
233 126 : CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
234 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, &
235 : -1.0_dp, dcdr_env%mo_coeff(ispin), dcdr_env%chc(ispin), &
236 126 : 0.0_dp, buf)
237 :
238 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, buf, dcdr_env%op_dR(ispin), &
239 126 : nmo, alpha=1.0_dp, beta=1.0_dp)
240 126 : CALL cp_fm_release(buf)
241 :
242 : ! SL multiply by -1 for response solver (H-S<H> C + dR_coupled= - (op_dR)
243 360 : CALL cp_fm_scale(-1.0_dp, dcdr_env%op_dR(ispin))
244 : END DO
245 :
246 108 : CALL dbcsr_deallocate_matrix_set(opdr_sym)
247 :
248 108 : CALL timestop(handle)
249 108 : END SUBROUTINE dcdr_build_op_dR
250 :
251 : ! **************************************************************************************************
252 : !> \brief Get the dC/dR by solving the Sternheimer equation, using the op_dR matrix
253 : !> \param dcdr_env ...
254 : !> \param p_env ...
255 : !> \param qs_env ...
256 : !> \authors SL, ED
257 : ! **************************************************************************************************
258 108 : SUBROUTINE dcdr_response_dR(dcdr_env, p_env, qs_env)
259 :
260 : TYPE(dcdr_env_type) :: dcdr_env
261 : TYPE(qs_p_env_type) :: p_env
262 : TYPE(qs_environment_type), POINTER :: qs_env
263 :
264 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_response_dR'
265 :
266 : INTEGER :: handle, ispin, output_unit
267 : LOGICAL :: should_stop
268 108 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h1_psi0, psi0_order, psi1
269 : TYPE(cp_fm_type), POINTER :: mo_coeff
270 : TYPE(cp_logger_type), POINTER :: logger
271 : TYPE(linres_control_type), POINTER :: linres_control
272 108 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
273 : TYPE(section_vals_type), POINTER :: lr_section
274 :
275 108 : CALL timeset(routineN, handle)
276 108 : NULLIFY (linres_control, lr_section, logger)
277 :
278 : CALL get_qs_env(qs_env=qs_env, &
279 : linres_control=linres_control, &
280 108 : mos=mos)
281 :
282 108 : logger => cp_get_default_logger()
283 108 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
284 :
285 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
286 108 : extension=".linresLog")
287 108 : IF (output_unit > 0) THEN
288 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
289 54 : "*** Self consistent optimization of the response wavefunction ***"
290 : END IF
291 :
292 : ! allocate the vectors
293 450 : ALLOCATE (psi0_order(dcdr_env%nspins))
294 450 : ALLOCATE (psi1(dcdr_env%nspins))
295 450 : ALLOCATE (h1_psi0(dcdr_env%nspins))
296 :
297 234 : DO ispin = 1, dcdr_env%nspins
298 126 : CALL cp_fm_create(psi1(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
299 126 : CALL cp_fm_create(h1_psi0(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
300 126 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
301 234 : psi0_order(ispin) = mo_coeff
302 : END DO
303 :
304 : ! Restart
305 108 : IF (linres_control%linres_restart) THEN
306 18 : CALL dcdr_read_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
307 : ELSE
308 198 : DO ispin = 1, dcdr_env%nspins
309 198 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
310 : END DO
311 : END IF
312 :
313 108 : IF (output_unit > 0) THEN
314 : WRITE (output_unit, "(T10,A,I4,A)") &
315 54 : "Response to the perturbation operator referring to atom ", dcdr_env%lambda, &
316 108 : " displaced in "//ACHAR(dcdr_env%beta + 119)
317 : END IF
318 234 : DO ispin = 1, dcdr_env%nspins
319 126 : CALL cp_fm_set_all(dcdr_env%dCR(ispin), 0.0_dp)
320 234 : CALL cp_fm_to_fm(dcdr_env%op_dR(ispin), h1_psi0(ispin))
321 : END DO
322 :
323 108 : linres_control%lr_triplet = .FALSE. ! we do singlet response
324 108 : linres_control%do_kernel = .TRUE.
325 108 : linres_control%converged = .FALSE.
326 :
327 : ! Position perturbation to get dCR
328 : ! (H0-E0) psi1 = (H1-E1) psi0
329 : ! psi1 = the perturbed wavefunction
330 : ! h1_psi0 = (H1-E1-S1*\varepsilon)
331 : ! psi0_order = the unperturbed wavefunction
332 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, &
333 108 : output_unit, should_stop)
334 234 : DO ispin = 1, dcdr_env%nspins
335 234 : CALL cp_fm_to_fm(psi1(ispin), dcdr_env%dCR(ispin))
336 : END DO
337 :
338 : ! Write the new result to the restart file
339 108 : IF (linres_control%linres_restart) THEN
340 18 : CALL dcdr_write_restart(qs_env, lr_section, psi1, dcdr_env%lambda, dcdr_env%beta, "dCdR")
341 : END IF
342 :
343 : ! clean up
344 234 : DO ispin = 1, dcdr_env%nspins
345 126 : CALL cp_fm_release(psi1(ispin))
346 234 : CALL cp_fm_release(h1_psi0(ispin))
347 : END DO
348 108 : DEALLOCATE (psi1, h1_psi0, psi0_order)
349 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
350 108 : "PRINT%PROGRAM_RUN_INFO")
351 :
352 108 : CALL timestop(handle)
353 :
354 324 : END SUBROUTINE dcdr_response_dR
355 :
356 : ! **************************************************************************************************
357 : !> \brief Calculate atomic polar tensor
358 : !> \param qs_env ...
359 : !> \param dcdr_env ...
360 : !> \author Edward Ditler
361 : ! **************************************************************************************************
362 216 : SUBROUTINE apt_dR(qs_env, dcdr_env)
363 : TYPE(qs_environment_type), POINTER :: qs_env
364 : TYPE(dcdr_env_type) :: dcdr_env
365 :
366 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR'
367 :
368 : INTEGER :: alpha, handle, ikind, ispin, nao, nmo
369 : LOGICAL :: ghost
370 : REAL(dp) :: apt_basis_derivative, &
371 : apt_coeff_derivative, charge, f_spin
372 72 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
373 : TYPE(cp_fm_type) :: overlap1_MO, tmp_fm_like_mos
374 : TYPE(cp_fm_type), POINTER :: mo_coeff
375 72 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
376 72 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
377 :
378 72 : apt_basis_derivative = 0._dp
379 72 : apt_coeff_derivative = 0._dp
380 :
381 72 : CALL timeset(routineN, handle)
382 :
383 72 : NULLIFY (qs_kind_set, particle_set)
384 : CALL get_qs_env(qs_env=qs_env, &
385 : qs_kind_set=qs_kind_set, &
386 72 : particle_set=particle_set)
387 :
388 72 : nao = dcdr_env%nao
389 72 : apt_el => dcdr_env%apt_el_dcdr
390 72 : apt_nuc => dcdr_env%apt_nuc_dcdr
391 :
392 72 : f_spin = 2._dp/dcdr_env%nspins
393 :
394 162 : DO ispin = 1, dcdr_env%nspins
395 : ! Compute S^(1,R)_(ij)
396 90 : CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
397 90 : CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
398 90 : nmo = dcdr_env%nmo(ispin)
399 90 : mo_coeff => dcdr_env%mo_coeff(ispin)
400 90 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
401 90 : CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
402 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
403 90 : tmp_fm_like_mos, ncol=nmo)
404 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
405 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
406 90 : 0.0_dp, overlap1_MO)
407 :
408 : ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
409 : ! We get the negative of the coefficients out of the linres solver
410 : ! And apply the constant correction due to the overlap derivative.
411 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
412 : -0.5_dp, mo_coeff, overlap1_MO, &
413 90 : -1.0_dp, dcdr_env%dCR_prime(ispin))
414 90 : CALL cp_fm_release(overlap1_MO)
415 :
416 360 : DO alpha = 1, 3
417 : ! FIRST CONTRIBUTION: dCR * moments * mo
418 270 : CALL cp_fm_set_all(tmp_fm_like_mos, 0._dp)
419 270 : CALL dbcsr_desymmetrize(dcdr_env%matrix_s1(1)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
420 270 : CALL dbcsr_desymmetrize(dcdr_env%moments(alpha)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix)
421 : CALL dbcsr_add(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%matrix_nosym_temp(2)%matrix, &
422 270 : -dcdr_env%ref_point(alpha), 1._dp)
423 :
424 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_nosym_temp(1)%matrix, dcdr_env%dCR_prime(ispin), &
425 270 : tmp_fm_like_mos, ncol=nmo)
426 270 : CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_coeff_derivative)
427 :
428 270 : apt_coeff_derivative = (-2._dp)*f_spin*apt_coeff_derivative
429 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
430 270 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
431 :
432 : ! SECOND CONTRIBUTION: We assemble all combinations of r_i, d(chi)/d(idir)
433 : ! difdip contains derivatives with respect to atom dcdr_env%lambda
434 : ! difdip(alpha, beta): < a | r_alpha | db/dR_beta >
435 : ! Multiply by the MO coefficients
436 270 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
437 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, mo_coeff, &
438 270 : tmp_fm_like_mos, ncol=nmo)
439 270 : CALL cp_fm_trace(mo_coeff, tmp_fm_like_mos, apt_basis_derivative)
440 :
441 : ! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
442 270 : apt_basis_derivative = -f_spin*apt_basis_derivative
443 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) = &
444 630 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
445 : END DO ! alpha
446 :
447 252 : CALL cp_fm_release(tmp_fm_like_mos)
448 : END DO !ispin
449 :
450 : ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
451 72 : CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
452 72 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
453 72 : IF (.NOT. ghost) THEN
454 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
455 72 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
456 : END IF
457 :
458 : ! And deallocate all the things!
459 72 : CALL cp_fm_release(tmp_fm_like_mos)
460 72 : CALL cp_fm_release(overlap1_MO)
461 :
462 72 : CALL timestop(handle)
463 72 : END SUBROUTINE apt_dR
464 :
465 : ! **************************************************************************************************
466 : !> \brief Calculate atomic polar tensor using the localized dipole operator
467 : !> \param qs_env ...
468 : !> \param dcdr_env ...
469 : !> \author Edward Ditler
470 : ! **************************************************************************************************
471 36 : SUBROUTINE apt_dR_localization(qs_env, dcdr_env)
472 : TYPE(qs_environment_type), POINTER :: qs_env
473 : TYPE(dcdr_env_type) :: dcdr_env
474 :
475 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apt_dR_localization'
476 :
477 : INTEGER :: alpha, handle, i, icenter, ikind, ispin, &
478 : map_atom, map_molecule, &
479 : max_nbr_center, nao, natom, nmo, &
480 : nsubset
481 : INTEGER, ALLOCATABLE, DIMENSION(:) :: mapping_atom_molecule
482 36 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mapping_wannier_atom
483 : LOGICAL :: ghost
484 : REAL(dp) :: apt_basis_derivative, &
485 : apt_coeff_derivative, charge, f_spin, &
486 : smallest_r, this_factor, tmp_aptcontr, &
487 : tmp_r
488 36 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: diagonal_elements
489 : REAL(dp), DIMENSION(3) :: distance, r_shifted
490 36 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
491 36 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
492 : TYPE(cell_type), POINTER :: cell
493 36 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
494 : TYPE(cp_fm_type), POINTER :: mo_coeff, overlap1_MO, tmp_fm, &
495 : tmp_fm_like_mos, tmp_fm_momo
496 36 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
497 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
498 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
499 :
500 36 : CALL timeset(routineN, handle)
501 :
502 36 : NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
503 :
504 : CALL get_qs_env(qs_env=qs_env, &
505 : qs_kind_set=qs_kind_set, &
506 : particle_set=particle_set, &
507 : molecule_set=molecule_set, &
508 36 : cell=cell)
509 :
510 36 : nsubset = SIZE(molecule_set)
511 36 : natom = SIZE(particle_set)
512 36 : apt_el => dcdr_env%apt_el_dcdr
513 36 : apt_nuc => dcdr_env%apt_nuc_dcdr
514 36 : apt_subset => dcdr_env%apt_el_dcdr_per_subset
515 36 : apt_center => dcdr_env%apt_el_dcdr_per_center
516 :
517 : ! Map wannier functions to atoms
518 36 : IF (dcdr_env%nspins == 1) THEN
519 36 : max_nbr_center = dcdr_env%nbr_center(1)
520 : ELSE
521 0 : max_nbr_center = MAX(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
522 : END IF
523 144 : ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
524 108 : ALLOCATE (mapping_atom_molecule(natom))
525 36 : centers_set => dcdr_env%centers_set
526 :
527 72 : DO ispin = 1, dcdr_env%nspins
528 180 : DO icenter = 1, dcdr_env%nbr_center(ispin)
529 : ! For every center we check which atom is closest
530 : CALL shift_wannier_into_cell(r=centers_set(ispin)%array(1:3, icenter), &
531 : cell=cell, &
532 144 : r_shifted=r_shifted)
533 :
534 144 : smallest_r = HUGE(0._dp)
535 612 : DO i = 1, natom
536 432 : distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
537 1728 : tmp_r = SUM(distance**2)
538 576 : IF (tmp_r < smallest_r) THEN
539 216 : mapping_wannier_atom(icenter, ispin) = i
540 216 : smallest_r = tmp_r
541 : END IF
542 : END DO
543 : END DO
544 :
545 : ! Map atoms to molecules
546 36 : CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
547 72 : IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
548 20 : DO icenter = 1, dcdr_env%nbr_center(ispin)
549 16 : map_atom = mapping_wannier_atom(icenter, ispin)
550 20 : map_molecule = mapping_atom_molecule(map_atom)
551 : END DO
552 : END IF
553 : END DO !ispin
554 :
555 36 : nao = dcdr_env%nao
556 36 : f_spin = 2._dp/dcdr_env%nspins
557 :
558 72 : DO ispin = 1, dcdr_env%nspins
559 : ! Compute S^(1,R)_(ij)
560 :
561 36 : ALLOCATE (tmp_fm_like_mos)
562 36 : ALLOCATE (overlap1_MO)
563 36 : CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
564 36 : CALL cp_fm_create(overlap1_MO, dcdr_env%momo_fm_struct(ispin)%struct)
565 36 : nmo = dcdr_env%nmo(ispin)
566 36 : mo_coeff => dcdr_env%mo_coeff(ispin)
567 36 : CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
568 36 : CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
569 : CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
570 36 : tmp_fm_like_mos, ncol=nmo)
571 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
572 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
573 36 : 0.0_dp, overlap1_MO)
574 :
575 : ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
576 : ! We get the negative of the coefficients out of the linres solver
577 : ! And apply the constant correction due to the overlap derivative.
578 : CALL parallel_gemm("N", "N", nao, nmo, nmo, &
579 : -0.5_dp, mo_coeff, overlap1_MO, &
580 36 : -1.0_dp, dcdr_env%dCR_prime(ispin))
581 36 : CALL cp_fm_release(overlap1_MO)
582 :
583 108 : ALLOCATE (diagonal_elements(nmo))
584 :
585 : ! Allocate temporary matrices
586 36 : ALLOCATE (tmp_fm)
587 36 : ALLOCATE (tmp_fm_momo)
588 36 : CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
589 36 : CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
590 :
591 : ! FIRST CONTRIBUTION: dCR * moments * mo
592 36 : this_factor = -2._dp*f_spin
593 144 : DO alpha = 1, 3
594 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
595 432 : CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
596 : CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
597 432 : ref_point=centers_set(ispin)%array(1:3, icenter))
598 : CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
599 : mo_coeff=dcdr_env%dCR_prime(ispin), work=tmp_fm, nmo=nmo, &
600 : icenter=icenter, &
601 540 : res=tmp_fm_like_mos)
602 : END DO
603 :
604 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
605 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
606 108 : 0.0_dp, tmp_fm_momo)
607 108 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
608 :
609 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
610 432 : map_atom = mapping_wannier_atom(icenter, ispin)
611 432 : map_molecule = mapping_atom_molecule(map_atom)
612 432 : tmp_aptcontr = this_factor*diagonal_elements(icenter)
613 :
614 : apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
615 432 : = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
616 :
617 : apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
618 540 : = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
619 : END DO
620 :
621 540 : apt_coeff_derivative = this_factor*SUM(diagonal_elements)
622 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
623 144 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_coeff_derivative
624 : END DO
625 :
626 : ! SECOND CONTRIBUTION: We assemble all combinations of r_i, dphi/d(idir)
627 : ! build part with AOs differentiated with respect to nuclear coordinates
628 : ! difdip contains derivatives with respect to atom dcdr_env%lambda
629 : ! difdip(alpha, beta): < a | r_alpha | d b/dR_beta >
630 36 : this_factor = -f_spin
631 144 : DO alpha = 1, 3
632 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
633 : ! Build the AO matrix with the right wannier center as reference point
634 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(1, dcdr_env%beta)%matrix, 0._dp)
635 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(2, dcdr_env%beta)%matrix, 0._dp)
636 432 : CALL dbcsr_set(dcdr_env%matrix_difdip(3, dcdr_env%beta)%matrix, 0._dp)
637 : CALL dipole_deriv_ao(qs_env, dcdr_env%matrix_difdip, dcdr_env%delta_basis_function, &
638 432 : 1, centers_set(ispin)%array(1:3, icenter))
639 : CALL multiply_localization(ao_matrix=dcdr_env%matrix_difdip(alpha, dcdr_env%beta)%matrix, &
640 : mo_coeff=mo_coeff, work=tmp_fm, nmo=nmo, &
641 : icenter=icenter, &
642 540 : res=tmp_fm_like_mos)
643 : END DO ! icenter
644 :
645 : CALL parallel_gemm("T", "N", nmo, nmo, nao, &
646 : 1.0_dp, mo_coeff, tmp_fm_like_mos, &
647 108 : 0.0_dp, tmp_fm_momo)
648 108 : CALL cp_fm_get_diag(tmp_fm_momo, diagonal_elements)
649 :
650 540 : DO icenter = 1, dcdr_env%nbr_center(ispin)
651 432 : map_atom = mapping_wannier_atom(icenter, ispin)
652 432 : map_molecule = mapping_atom_molecule(map_atom)
653 432 : tmp_aptcontr = this_factor*diagonal_elements(icenter)
654 :
655 : apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) &
656 432 : = apt_subset(dcdr_env%beta, alpha, dcdr_env%lambda, map_molecule) + tmp_aptcontr
657 :
658 : apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) &
659 540 : = apt_center(dcdr_env%beta, alpha, dcdr_env%lambda, icenter) + tmp_aptcontr
660 : END DO
661 :
662 : ! The negative sign is due to dipole_deriv_ao computing the derivatives with respect to nuclear coordinates.
663 540 : apt_basis_derivative = this_factor*SUM(diagonal_elements)
664 :
665 : apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) &
666 144 : = apt_el(dcdr_env%beta, alpha, dcdr_env%lambda) + apt_basis_derivative
667 :
668 : END DO ! alpha
669 36 : DEALLOCATE (diagonal_elements)
670 :
671 36 : CALL cp_fm_release(tmp_fm)
672 36 : CALL cp_fm_release(tmp_fm_like_mos)
673 36 : CALL cp_fm_release(tmp_fm_momo)
674 36 : DEALLOCATE (overlap1_MO)
675 36 : DEALLOCATE (tmp_fm)
676 36 : DEALLOCATE (tmp_fm_like_mos)
677 144 : DEALLOCATE (tmp_fm_momo)
678 : END DO !ispin
679 :
680 : ! Finally the nuclear contribution: nuclear charge * Kronecker_delta_{dcdr_env%beta,i}
681 36 : CALL get_atomic_kind(particle_set(dcdr_env%lambda)%atomic_kind, kind_number=ikind)
682 36 : CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost)
683 36 : IF (.NOT. ghost) THEN
684 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) = &
685 36 : apt_nuc(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda) + charge
686 :
687 36 : map_molecule = mapping_atom_molecule(dcdr_env%lambda)
688 : apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) &
689 36 : = apt_subset(dcdr_env%beta, dcdr_env%beta, dcdr_env%lambda, map_molecule) + charge
690 : END IF
691 :
692 : ! And deallocate all the things!
693 :
694 36 : CALL timestop(handle)
695 108 : END SUBROUTINE apt_dR_localization
696 :
697 : END MODULE qs_dcdr
|
