Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate the derivatives of the MO coefficients wrt nuclear coordinates
10 : !> \author Sandra Luber, Edward Ditler
11 : ! **************************************************************************************************
12 :
13 : MODULE qs_dcdr_utils
14 : !#include "./common/cp_common_uses.f90"
15 : USE cell_types, ONLY: cell_type,&
16 : get_cell
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
19 : dbcsr_create,&
20 : dbcsr_init_p,&
21 : dbcsr_p_type,&
22 : dbcsr_set,&
23 : dbcsr_type,&
24 : dbcsr_type_no_symmetry,&
25 : dbcsr_type_symmetric
26 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
27 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
28 : dbcsr_allocate_matrix_set,&
29 : dbcsr_deallocate_matrix_set
30 : USE cp_files, ONLY: close_file,&
31 : open_file
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release
34 : USE cp_fm_types, ONLY: cp_fm_create,&
35 : cp_fm_get_info,&
36 : cp_fm_get_submatrix,&
37 : cp_fm_release,&
38 : cp_fm_set_all,&
39 : cp_fm_set_submatrix,&
40 : cp_fm_to_fm,&
41 : cp_fm_type
42 : USE cp_log_handling, ONLY: cp_get_default_logger,&
43 : cp_logger_get_default_io_unit,&
44 : cp_logger_type,&
45 : cp_to_string
46 : USE cp_output_handling, ONLY: cp_p_file,&
47 : cp_print_key_finished_output,&
48 : cp_print_key_generate_filename,&
49 : cp_print_key_should_output,&
50 : cp_print_key_unit_nr
51 : USE cp_result_methods, ONLY: get_results
52 : USE cp_result_types, ONLY: cp_result_type
53 : USE input_constants, ONLY: current_orb_center_wannier,&
54 : use_mom_ref_user
55 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
56 : section_vals_type,&
57 : section_vals_val_get
58 : USE kinds, ONLY: default_path_length,&
59 : default_string_length,&
60 : dp
61 : USE memory_utilities, ONLY: reallocate
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE molecule_types, ONLY: molecule_type
64 : USE moments_utils, ONLY: get_reference_point
65 : USE parallel_gemm_api, ONLY: parallel_gemm
66 : USE particle_types, ONLY: particle_type
67 : USE qs_environment_types, ONLY: get_qs_env,&
68 : qs_environment_type
69 : USE qs_kinetic, ONLY: build_kinetic_matrix
70 : USE qs_ks_types, ONLY: qs_ks_env_type
71 : USE qs_linres_types, ONLY: dcdr_env_type,&
72 : linres_control_type
73 : USE qs_loc_types, ONLY: get_qs_loc_env,&
74 : localized_wfn_control_type,&
75 : qs_loc_env_type
76 : USE qs_mo_types, ONLY: get_mo_set,&
77 : mo_set_type
78 : USE qs_moments, ONLY: build_local_moment_matrix
79 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
80 : USE qs_overlap, ONLY: build_overlap_matrix
81 : USE string_utilities, ONLY: xstring
82 : #include "./base/base_uses.f90"
83 :
84 : IMPLICIT NONE
85 :
86 : PRIVATE
87 : PUBLIC :: dcdr_env_cleanup, dcdr_env_init, dcdr_print, &
88 : get_loc_setting, shift_wannier_into_cell, &
89 : dcdr_write_restart, dcdr_read_restart, &
90 : multiply_localization
91 :
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_utils'
93 :
94 : REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = RESHAPE((/ &
95 : 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
96 : 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
97 : 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
98 : (/3, 3, 3/))
99 :
100 : CONTAINS
101 :
102 : ! **************************************************************************************************
103 : !> \brief Multiply (ao_matrix @ mo_coeff) and store the column icenter in res
104 : !> \param ao_matrix ...
105 : !> \param mo_coeff ...
106 : !> \param work Working space
107 : !> \param nmo ...
108 : !> \param icenter ...
109 : !> \param res ...
110 : !> \author Edward Ditler
111 : ! **************************************************************************************************
112 864 : SUBROUTINE multiply_localization(ao_matrix, mo_coeff, work, nmo, icenter, res)
113 : TYPE(dbcsr_type), INTENT(IN), POINTER :: ao_matrix
114 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
115 : TYPE(cp_fm_type), INTENT(INOUT) :: work
116 : INTEGER, INTENT(IN) :: nmo, icenter
117 : TYPE(cp_fm_type), INTENT(IN) :: res
118 :
119 : CHARACTER(LEN=*), PARAMETER :: routineN = 'multiply_localization'
120 :
121 : INTEGER :: handle
122 :
123 864 : CALL timeset(routineN, handle)
124 :
125 : ! Multiply by the MO coefficients
126 864 : CALL cp_dbcsr_sm_fm_multiply(ao_matrix, mo_coeff, work, ncol=nmo)
127 :
128 : ! Only keep the icenter-th column
129 864 : CALL cp_fm_to_fm(work, res, 1, icenter, icenter)
130 :
131 : ! Reset the matrices
132 864 : CALL cp_fm_set_all(work, 0.0_dp)
133 :
134 864 : CALL timestop(handle)
135 864 : END SUBROUTINE multiply_localization
136 :
137 : ! **************************************************************************************************
138 : !> \brief Copied from linres_read_restart
139 : !> \param qs_env ...
140 : !> \param linres_section ...
141 : !> \param vec ...
142 : !> \param lambda ...
143 : !> \param beta ...
144 : !> \param tag ...
145 : !> \note Adapted from linres_read_restart (ED)
146 : !> Would be nice not to crash but to start from zero if the present file doesn't match.
147 : ! **************************************************************************************************
148 18 : SUBROUTINE dcdr_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
149 : TYPE(qs_environment_type), POINTER :: qs_env
150 : TYPE(section_vals_type), POINTER :: linres_section
151 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
152 : INTEGER, INTENT(IN) :: lambda, beta
153 : CHARACTER(LEN=*) :: tag
154 :
155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_read_restart'
156 :
157 : CHARACTER(LEN=default_path_length) :: filename
158 : CHARACTER(LEN=default_string_length) :: my_middle
159 : INTEGER :: beta_tmp, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, lambda_tmp, &
160 : max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
161 : LOGICAL :: file_exists
162 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
163 : TYPE(cp_fm_type), POINTER :: mo_coeff
164 : TYPE(cp_logger_type), POINTER :: logger
165 18 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
166 : TYPE(mp_para_env_type), POINTER :: para_env
167 : TYPE(section_vals_type), POINTER :: print_key
168 :
169 18 : file_exists = .FALSE.
170 :
171 18 : CALL timeset(routineN, handle)
172 :
173 18 : NULLIFY (mos, para_env, logger, print_key, vecbuffer)
174 18 : logger => cp_get_default_logger()
175 :
176 : iounit = cp_print_key_unit_nr(logger, linres_section, &
177 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
178 :
179 : CALL get_qs_env(qs_env=qs_env, &
180 : para_env=para_env, &
181 18 : mos=mos)
182 :
183 18 : nspins = SIZE(mos)
184 :
185 18 : rst_unit = -1
186 18 : IF (para_env%is_source()) THEN
187 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
188 9 : n_rep_val=n_rep_val)
189 :
190 9 : CALL XSTRING(tag, ia, ie)
191 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
192 9 : //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
193 :
194 9 : IF (n_rep_val > 0) THEN
195 0 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
196 0 : CALL xstring(filename, ia, ie)
197 0 : filename = filename(ia:ie)//TRIM(my_middle)//".lr"
198 : ELSE
199 : ! try to read from the filename that is generated automatically from the printkey
200 9 : print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
201 : filename = cp_print_key_generate_filename(logger, print_key, &
202 9 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
203 : END IF
204 9 : INQUIRE (FILE=filename, exist=file_exists)
205 : !
206 : ! open file
207 9 : IF (file_exists) THEN
208 : CALL open_file(file_name=TRIM(filename), &
209 : file_action="READ", &
210 : file_form="UNFORMATTED", &
211 : file_position="REWIND", &
212 : file_status="OLD", &
213 0 : unit_number=rst_unit)
214 :
215 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
216 0 : "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
217 : ELSE
218 9 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
219 9 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
220 : END IF
221 : END IF
222 :
223 18 : CALL para_env%bcast(file_exists)
224 :
225 18 : IF (file_exists) THEN
226 :
227 0 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
228 0 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
229 :
230 0 : ALLOCATE (vecbuffer(nao, max_block))
231 : !
232 : ! read headers
233 0 : IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
234 0 : CALL para_env%bcast(iostat)
235 0 : CALL para_env%bcast(beta_tmp)
236 0 : CALL para_env%bcast(lambda_tmp)
237 0 : CALL para_env%bcast(nspins_tmp)
238 0 : CALL para_env%bcast(nao_tmp)
239 :
240 : ! check that the number nao, nmo and nspins are
241 : ! the same as in the current mos
242 0 : IF (nspins_tmp .NE. nspins) THEN
243 0 : CPABORT("nspins not consistent")
244 : END IF
245 0 : IF (nao_tmp .NE. nao) CPABORT("nao not consistent")
246 : ! check that it's the right file
247 : ! the same as in the current mos
248 0 : IF (lambda_tmp .NE. lambda) CPABORT("lambda not consistent")
249 0 : IF (beta_tmp .NE. beta) CPABORT("beta not consistent")
250 : !
251 0 : DO ispin = 1, nspins
252 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
253 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
254 : !
255 0 : IF (rst_unit > 0) READ (rst_unit) nmo_tmp
256 0 : CALL para_env%bcast(nmo_tmp)
257 0 : IF (nmo_tmp .NE. nmo) CPABORT("nmo not consistent")
258 : !
259 : ! read the response
260 0 : DO i = 1, nmo, MAX(max_block, 1)
261 0 : i_block = MIN(max_block, nmo - i + 1)
262 0 : DO j = 1, i_block
263 0 : IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
264 : END DO
265 0 : CALL para_env%bcast(vecbuffer)
266 0 : CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
267 : END DO
268 : END DO
269 :
270 0 : IF (iostat /= 0) THEN
271 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
272 0 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
273 : END IF
274 :
275 0 : DEALLOCATE (vecbuffer)
276 :
277 : END IF
278 :
279 18 : IF (para_env%is_source()) THEN
280 9 : IF (file_exists) CALL close_file(unit_number=rst_unit)
281 : END IF
282 :
283 18 : CALL timestop(handle)
284 :
285 18 : END SUBROUTINE dcdr_read_restart
286 :
287 : ! **************************************************************************************************
288 : !> \brief Copied from linres_write_restart
289 : !> \param qs_env ...
290 : !> \param linres_section ...
291 : !> \param vec ...
292 : !> \param lambda ...
293 : !> \param beta ...
294 : !> \param tag ...
295 : !> \note Adapted from linres_read_restart (ED)
296 : !> Would be nice not to crash but to start from zero if the present file doesn't match.
297 : ! **************************************************************************************************
298 18 : SUBROUTINE dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
299 : TYPE(qs_environment_type), POINTER :: qs_env
300 : TYPE(section_vals_type), POINTER :: linres_section
301 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
302 : INTEGER, INTENT(IN) :: lambda, beta
303 : CHARACTER(LEN=*) :: tag
304 :
305 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_write_restart'
306 :
307 : CHARACTER(LEN=default_path_length) :: filename
308 : CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status
309 : INTEGER :: handle, i, i_block, ia, ie, iounit, &
310 : ispin, j, max_block, nao, nmo, nspins, &
311 : rst_unit
312 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
313 : TYPE(cp_fm_type), POINTER :: mo_coeff
314 : TYPE(cp_logger_type), POINTER :: logger
315 18 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
316 : TYPE(mp_para_env_type), POINTER :: para_env
317 : TYPE(section_vals_type), POINTER :: print_key
318 :
319 18 : NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
320 :
321 18 : CALL timeset(routineN, handle)
322 :
323 18 : logger => cp_get_default_logger()
324 :
325 18 : IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
326 : used_print_key=print_key), &
327 : cp_p_file)) THEN
328 :
329 : iounit = cp_print_key_unit_nr(logger, linres_section, &
330 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
331 :
332 : CALL get_qs_env(qs_env=qs_env, &
333 : mos=mos, &
334 18 : para_env=para_env)
335 :
336 18 : nspins = SIZE(mos)
337 :
338 18 : my_status = "REPLACE"
339 18 : my_pos = "REWIND"
340 18 : CALL XSTRING(tag, ia, ie)
341 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
342 18 : //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
343 : rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
344 : extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), &
345 18 : file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
346 :
347 : filename = cp_print_key_generate_filename(logger, print_key, &
348 18 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
349 :
350 18 : IF (iounit > 0) THEN
351 : WRITE (UNIT=iounit, FMT="(T2,A)") &
352 9 : "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">"
353 : END IF
354 :
355 : !
356 : ! write data to file
357 : ! use the scalapack block size as a default for buffering columns
358 18 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
359 18 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
360 72 : ALLOCATE (vecbuffer(nao, max_block))
361 :
362 18 : IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
363 :
364 36 : DO ispin = 1, nspins
365 18 : CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
366 :
367 18 : IF (rst_unit > 0) WRITE (rst_unit) nmo
368 :
369 72 : DO i = 1, nmo, MAX(max_block, 1)
370 18 : i_block = MIN(max_block, nmo - i + 1)
371 18 : CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
372 : ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
373 : ! to old ones, and in cases where max_block is different between runs, as might happen during
374 : ! restarts with a different number of CPUs
375 108 : DO j = 1, i_block
376 90 : IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
377 : END DO
378 : END DO
379 : END DO
380 :
381 18 : DEALLOCATE (vecbuffer)
382 :
383 : CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
384 36 : "PRINT%RESTART")
385 : END IF
386 :
387 18 : CALL timestop(handle)
388 :
389 18 : END SUBROUTINE dcdr_write_restart
390 :
391 : ! **************************************************************************************************
392 : !> \brief Print the APT and sum rules
393 : !> \param dcdr_env ...
394 : !> \param qs_env ...
395 : !> \author Edward Ditler
396 : ! **************************************************************************************************
397 22 : SUBROUTINE dcdr_print(dcdr_env, qs_env)
398 : TYPE(dcdr_env_type) :: dcdr_env
399 : TYPE(qs_environment_type), POINTER :: qs_env
400 :
401 : CHARACTER(LEN=default_string_length) :: description
402 : INTEGER :: alpha, beta, delta, gamma, i, k, l, &
403 : lambda, natom, nsubset, output_unit
404 22 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
405 22 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center_dcdr, apt_subset_dcdr
406 : REAL(kind=dp), DIMENSION(3, 3) :: sum_rule_0, sum_rule_1, sum_rule_2
407 : TYPE(cp_logger_type), POINTER :: logger
408 : TYPE(cp_result_type), POINTER :: results
409 22 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
410 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
411 : TYPE(section_vals_type), POINTER :: dcdr_section
412 :
413 22 : NULLIFY (logger)
414 :
415 44 : logger => cp_get_default_logger()
416 22 : output_unit = cp_logger_get_default_io_unit(logger)
417 :
418 22 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
419 :
420 22 : NULLIFY (particle_set)
421 22 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
422 22 : natom = SIZE(particle_set)
423 22 : nsubset = SIZE(molecule_set)
424 :
425 22 : apt_el_dcdr => dcdr_env%apt_el_dcdr
426 22 : apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr
427 22 : apt_total_dcdr => dcdr_env%apt_total_dcdr
428 22 : apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset
429 22 : apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center
430 :
431 22 : IF (dcdr_env%localized_psi0) THEN
432 4 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the final apt matrix per atom per subset'
433 16 : DO k = 1, natom
434 52 : DO l = 1, nsubset
435 36 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, *) 'APT | Subset', l
436 156 : DO i = 1, 3
437 108 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,F15.6,F15.6,F15.6)") &
438 90 : 'APT | apt_subset ', i, apt_subset_dcdr(i, :, k, l)
439 : END DO
440 : END DO
441 : END DO
442 : END IF
443 :
444 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") &
445 11 : 'APT | Write the final apt matrix per atom (Position perturbation)'
446 88 : DO l = 1, natom
447 66 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") &
448 33 : 'APT | Atom', l, ' - GAPT ', &
449 : (apt_total_dcdr(1, 1, l) &
450 : + apt_total_dcdr(2, 2, l) &
451 66 : + apt_total_dcdr(3, 3, l))/3._dp
452 286 : DO i = 1, 3
453 264 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
454 : END DO
455 : END DO
456 :
457 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix'
458 88 : DO i = 1, 3
459 66 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
460 451 : "(A,F15.6,F15.6,F15.6)") "APT | ", SUM(apt_total_dcdr(i, :, :), dim=2)
461 : END DO
462 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | End Write the final apt matrix'
463 :
464 : ! Get the dipole
465 22 : CALL get_qs_env(qs_env, results=results)
466 22 : description = "[DIPOLE]"
467 22 : CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3))
468 :
469 : ! Sum rules [for all alpha, beta]
470 22 : sum_rule_0 = 0._dp
471 22 : sum_rule_1 = 0._dp
472 22 : sum_rule_2 = 0._dp
473 :
474 88 : DO alpha = 1, 3
475 286 : DO beta = 1, 3
476 : ! 0: sum_lambda apt(alpha, beta, lambda)
477 792 : DO lambda = 1, natom
478 : sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
479 792 : + apt_total_dcdr(alpha, beta, lambda)
480 : END DO
481 :
482 : ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
483 792 : DO gamma = 1, 3
484 : sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
485 792 : + Levi_Civita(alpha, beta, gamma)*dcdr_env%dipole_pos(gamma)
486 : END DO
487 :
488 : ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
489 858 : DO lambda = 1, natom
490 2574 : DO gamma = 1, 3
491 7722 : DO delta = 1, 3
492 : sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
493 : + Levi_Civita(beta, gamma, delta) &
494 : *particle_set(lambda)%r(gamma) &
495 7128 : *apt_total_dcdr(delta, alpha, lambda)
496 : END DO
497 : END DO
498 : END DO
499 :
500 : END DO ! beta
501 : END DO ! alpha
502 :
503 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules"
504 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") &
505 11 : "APT |", "Total APT", "Dipole", "R * APT"
506 88 : DO alpha = 1, 3
507 286 : DO beta = 1, 3
508 198 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
509 : "(A,I3,I3,F15.6,F15.6,F15.6)") &
510 99 : "APT | ", &
511 99 : alpha, beta, &
512 99 : sum_rule_0(alpha, beta), &
513 99 : sum_rule_1(alpha, beta), &
514 264 : sum_rule_2(alpha, beta)
515 : END DO
516 : END DO
517 :
518 22 : END SUBROUTINE dcdr_print
519 :
520 : ! **************************************************************************************************
521 : !> \brief ...
522 : !> \param r ...
523 : !> \param cell ...
524 : !> \param r_shifted ...
525 : ! **************************************************************************************************
526 144 : SUBROUTINE shift_wannier_into_cell(r, cell, r_shifted)
527 : REAL(dp), DIMENSION(3), INTENT(in) :: r
528 : TYPE(cell_type), INTENT(in), POINTER :: cell
529 : REAL(dp), DIMENSION(3), INTENT(out) :: r_shifted
530 :
531 : INTEGER :: i
532 : REAL(kind=dp), DIMENSION(3) :: abc
533 :
534 : ! Only orthorombic cell for now
535 144 : CALL get_cell(cell, abc=abc)
536 :
537 576 : DO i = 1, 3
538 576 : IF (r(i) < 0._dp) THEN
539 234 : r_shifted(i) = r(i) + abc(i)
540 198 : ELSE IF (r(i) > abc(i)) THEN
541 0 : r_shifted(i) = r(i) - abc(i)
542 : ELSE
543 198 : r_shifted(i) = r(i)
544 : END IF
545 : END DO
546 144 : END SUBROUTINE shift_wannier_into_cell
547 :
548 : ! **************************************************************************************************
549 : !> \brief ...
550 : !> \param dcdr_env ...
551 : !> \param qs_env ...
552 : ! **************************************************************************************************
553 4 : SUBROUTINE get_loc_setting(dcdr_env, qs_env)
554 : TYPE(dcdr_env_type) :: dcdr_env
555 : TYPE(qs_environment_type), POINTER :: qs_env
556 :
557 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_loc_setting'
558 :
559 : INTEGER :: handle, is, ispin, istate, max_states, &
560 : nmo, nmoloc, nstate, nstate_list(2)
561 4 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
562 4 : REAL(dp), DIMENSION(:, :), POINTER :: center_array
563 : TYPE(linres_control_type), POINTER :: linres_control
564 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
565 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
566 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
567 : TYPE(section_vals_type), POINTER :: dcdr_section
568 :
569 4 : CALL timeset(routineN, handle)
570 :
571 : CALL get_qs_env(qs_env=qs_env, &
572 : linres_control=linres_control, &
573 4 : mos=mos)
574 :
575 : ! Some checks
576 4 : max_states = 0
577 4 : CALL get_mo_set(mo_set=mos(1), nmo=nmo)
578 4 : max_states = MAX(max_states, nmo)
579 :
580 : ! check that the number of localized states is equal to the number of states
581 4 : nmoloc = SIZE(linres_control%qs_loc_env%localized_wfn_control%centers_set(1)%array, 2)
582 4 : IF (nmoloc .NE. nmo) THEN
583 0 : CPABORT("The number of localized functions is not equal to the number of states.")
584 : END IF
585 :
586 : ! which center for the orbitals shall we use
587 4 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
588 4 : CALL section_vals_val_get(dcdr_section, "ORBITAL_CENTER", i_val=dcdr_env%orb_center)
589 8 : SELECT CASE (dcdr_env%orb_center)
590 : CASE (current_orb_center_wannier)
591 4 : dcdr_env%orb_center_name = "WANNIER"
592 : CASE DEFAULT
593 4 : CPABORT(" ")
594 : END SELECT
595 :
596 4 : qs_loc_env => linres_control%qs_loc_env
597 4 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
598 :
599 16 : ALLOCATE (dcdr_env%centers_set(dcdr_env%nspins))
600 12 : ALLOCATE (dcdr_env%center_list(dcdr_env%nspins))
601 16 : ALLOCATE (state_list(max_states, dcdr_env%nspins))
602 24 : state_list(:, :) = HUGE(0)
603 12 : nstate_list(:) = HUGE(0)
604 :
605 : ! Build the state_list
606 8 : DO ispin = 1, dcdr_env%nspins
607 4 : center_array => localized_wfn_control%centers_set(ispin)%array
608 4 : nstate = 0
609 20 : DO istate = 1, SIZE(center_array, 2)
610 16 : nstate = nstate + 1
611 20 : state_list(nstate, ispin) = istate
612 : END DO
613 : nstate_list(ispin) = nstate
614 :
615 : ! clustering the states
616 4 : nstate = nstate_list(ispin)
617 4 : dcdr_env%nstates(ispin) = nstate
618 :
619 12 : ALLOCATE (dcdr_env%center_list(ispin)%array(2, nstate + 1))
620 12 : ALLOCATE (dcdr_env%centers_set(ispin)%array(3, nstate))
621 64 : dcdr_env%center_list(ispin)%array(:, :) = HUGE(0)
622 68 : dcdr_env%centers_set(ispin)%array(:, :) = HUGE(0.0_dp)
623 :
624 4 : center_array => localized_wfn_control%centers_set(ispin)%array
625 :
626 : ! point to the psi0 centers
627 4 : SELECT CASE (dcdr_env%orb_center)
628 : CASE (current_orb_center_wannier)
629 : ! use the wannier center as -center-
630 4 : dcdr_env%nbr_center(ispin) = nstate
631 20 : DO is = 1, nstate
632 16 : istate = state_list(is, 1)
633 64 : dcdr_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
634 16 : dcdr_env%center_list(ispin)%array(1, is) = is
635 20 : dcdr_env%center_list(ispin)%array(2, is) = istate
636 : END DO
637 4 : dcdr_env%center_list(ispin)%array(1, nstate + 1) = nstate + 1
638 :
639 : CASE DEFAULT
640 4 : CPABORT("Unknown orbital center...")
641 : END SELECT
642 : END DO
643 :
644 4 : CALL timestop(handle)
645 12 : END SUBROUTINE get_loc_setting
646 :
647 : ! **************************************************************************************************
648 : !> \brief Initialize the dcdr environment
649 : !> \param dcdr_env ...
650 : !> \param qs_env ...
651 : ! **************************************************************************************************
652 24 : SUBROUTINE dcdr_env_init(dcdr_env, qs_env)
653 : TYPE(dcdr_env_type) :: dcdr_env
654 : TYPE(qs_environment_type), POINTER :: qs_env
655 :
656 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_env_init'
657 :
658 : INTEGER :: handle, homo, i, isize, ispin, j, jg, &
659 : n_rep, nao, natom, nmo, nspins, &
660 : nsubset, output_unit, reference, &
661 : unit_number
662 24 : INTEGER, DIMENSION(:), POINTER :: tmplist
663 : LOGICAL :: explicit
664 24 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
665 : TYPE(cp_fm_type) :: buf
666 : TYPE(cp_fm_type), POINTER :: mo_coeff
667 : TYPE(cp_logger_type), POINTER :: logger
668 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
669 : TYPE(dft_control_type), POINTER :: dft_control
670 24 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
671 24 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
672 : TYPE(mp_para_env_type), POINTER :: para_env
673 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
674 24 : POINTER :: sab_all, sab_orb
675 24 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
676 : TYPE(qs_ks_env_type), POINTER :: ks_env
677 : TYPE(section_vals_type), POINTER :: dcdr_section, loc_section, lr_section
678 :
679 24 : CALL timeset(routineN, handle)
680 :
681 : ! Set up the logger
682 24 : NULLIFY (logger, loc_section, dcdr_section, lr_section)
683 24 : logger => cp_get_default_logger()
684 24 : loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE")
685 24 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
686 : dcdr_env%output_unit = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", &
687 : extension=".data", middle_name="dcdr", log_filename=.FALSE., &
688 24 : file_position="REWIND", file_status="REPLACE")
689 :
690 24 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
691 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
692 24 : extension=".linresLog")
693 24 : unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
694 :
695 24 : IF (output_unit > 0) THEN
696 12 : WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***"
697 : END IF
698 :
699 24 : NULLIFY (ks_env, dft_control, sab_orb, sab_all, particle_set, molecule_set, matrix_s, matrix_ks, mos, para_env)
700 : CALL get_qs_env(qs_env=qs_env, &
701 : ks_env=ks_env, &
702 : dft_control=dft_control, &
703 : sab_orb=sab_orb, &
704 : sab_all=sab_all, &
705 : particle_set=particle_set, &
706 : molecule_set=molecule_set, &
707 : matrix_s=matrix_s, &
708 : matrix_ks=matrix_ks, &
709 : mos=mos, &
710 24 : para_env=para_env)
711 :
712 24 : natom = SIZE(particle_set)
713 24 : nsubset = SIZE(molecule_set)
714 24 : nspins = dft_control%nspins
715 24 : dcdr_env%nspins = dft_control%nspins
716 :
717 24 : NULLIFY (dcdr_env%matrix_s)
718 : CALL build_overlap_matrix(ks_env, matrix_s=dcdr_env%matrix_s, &
719 : matrix_name="OVERLAP MATRIX", &
720 : nderivative=1, &
721 : basis_type_a="ORB", &
722 : basis_type_b="ORB", &
723 24 : sab_nl=sab_orb)
724 :
725 24 : NULLIFY (dcdr_env%matrix_t)
726 : CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, &
727 : matrix_name="KINETIC ENERGY MATRIX", &
728 : basis_type="ORB", &
729 : sab_nl=sab_orb, nderivative=1, &
730 24 : eps_filter=dft_control%qs_control%eps_filter_matrix)
731 :
732 : ! Get inputs
733 24 : CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
734 24 : CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
735 24 : CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
736 24 : CALL section_vals_val_get(dcdr_section, "Z_MATRIX_METHOD", l_val=dcdr_env%z_matrix_method)
737 :
738 96 : dcdr_env%ref_point = 0._dp
739 :
740 : ! List of atoms
741 24 : NULLIFY (tmplist)
742 24 : isize = 0
743 24 : CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep)
744 24 : IF (n_rep == 0) THEN
745 72 : ALLOCATE (dcdr_env%list_of_atoms(natom))
746 96 : DO jg = 1, natom
747 96 : dcdr_env%list_of_atoms(jg) = jg
748 : END DO
749 : ELSE
750 0 : DO jg = 1, n_rep
751 0 : ALLOCATE (dcdr_env%list_of_atoms(isize))
752 0 : CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist)
753 0 : CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist))
754 0 : dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
755 0 : isize = SIZE(dcdr_env%list_of_atoms)
756 : END DO
757 : END IF
758 :
759 : ! Reference point
760 24 : IF (dcdr_env%localized_psi0) THEN
761 : ! Get the Wannier localized wave functions and centers
762 4 : CALL get_loc_setting(dcdr_env, qs_env)
763 : ELSE
764 : ! Get the reference point from the input
765 20 : CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference)
766 20 : CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit)
767 20 : IF (explicit) THEN
768 0 : CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point)
769 : ELSE
770 20 : IF (reference == use_mom_ref_user) &
771 0 : CPABORT("User-defined reference point should be given explicitly")
772 : END IF
773 :
774 : CALL get_reference_point(rpoint=dcdr_env%ref_point, qs_env=qs_env, &
775 : reference=reference, &
776 20 : ref_point=ref_point)
777 : END IF
778 :
779 : ! Helper matrix structs
780 : NULLIFY (dcdr_env%aoao_fm_struct, &
781 24 : dcdr_env%momo_fm_struct, &
782 24 : dcdr_env%likemos_fm_struct, &
783 24 : dcdr_env%homohomo_fm_struct)
784 : CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, &
785 24 : nao=nao, nmo=nmo, homo=homo)
786 : CALL cp_fm_struct_create(dcdr_env%aoao_fm_struct, nrow_global=nao, &
787 : ncol_global=nao, para_env=para_env, &
788 24 : context=mo_coeff%matrix_struct%context)
789 : CALL cp_fm_struct_create(dcdr_env%homohomo_fm_struct, nrow_global=homo, &
790 : ncol_global=homo, para_env=para_env, &
791 24 : context=mo_coeff%matrix_struct%context)
792 24 : dcdr_env%nao = nao
793 :
794 72 : ALLOCATE (dcdr_env%nmo(nspins))
795 100 : ALLOCATE (dcdr_env%momo_fm_struct(nspins))
796 76 : ALLOCATE (dcdr_env%likemos_fm_struct(nspins))
797 52 : DO ispin = 1, nspins
798 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
799 28 : nao=nao, nmo=nmo, homo=homo)
800 : CALL cp_fm_struct_create(dcdr_env%momo_fm_struct(ispin)%struct, nrow_global=nmo, &
801 : ncol_global=nmo, para_env=para_env, &
802 28 : context=mo_coeff%matrix_struct%context)
803 : CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct(ispin)%struct, &
804 28 : template_fmstruct=mo_coeff%matrix_struct)
805 80 : dcdr_env%nmo(ispin) = nmo
806 : END DO
807 :
808 : ! Fields of reals
809 72 : ALLOCATE (dcdr_env%deltaR(3, natom))
810 48 : ALLOCATE (dcdr_env%delta_basis_function(3, natom))
811 72 : ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom))
812 48 : ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom))
813 48 : ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom))
814 :
815 960 : dcdr_env%apt_el_dcdr = 0._dp
816 960 : dcdr_env%apt_nuc_dcdr = 0._dp
817 960 : dcdr_env%apt_total_dcdr = 0._dp
818 :
819 312 : dcdr_env%deltaR = 0.0_dp
820 312 : dcdr_env%delta_basis_function = 0._dp
821 :
822 : ! Localization
823 24 : IF (dcdr_env%localized_psi0) THEN
824 16 : ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1)))
825 16 : ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset))
826 12 : ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset))
827 644 : dcdr_env%apt_el_dcdr_per_center = 0._dp
828 484 : dcdr_env%apt_el_dcdr_per_subset = 0._dp
829 484 : dcdr_env%apt_subset = 0.0_dp
830 : END IF
831 :
832 : ! Full matrices
833 100 : ALLOCATE (dcdr_env%mo_coeff(nspins))
834 76 : ALLOCATE (dcdr_env%dCR(nspins))
835 76 : ALLOCATE (dcdr_env%dCR_prime(nspins))
836 76 : ALLOCATE (dcdr_env%chc(nspins))
837 76 : ALLOCATE (dcdr_env%op_dR(nspins))
838 :
839 52 : DO ispin = 1, nspins
840 28 : CALL cp_fm_create(dcdr_env%dCR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
841 28 : CALL cp_fm_create(dcdr_env%dCR_prime(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
842 28 : CALL cp_fm_create(dcdr_env%mo_coeff(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
843 28 : CALL cp_fm_create(dcdr_env%chc(ispin), dcdr_env%momo_fm_struct(ispin)%struct)
844 28 : CALL cp_fm_create(dcdr_env%op_dR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
845 :
846 28 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
847 52 : CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
848 : END DO
849 :
850 24 : IF (dcdr_env%z_matrix_method) THEN
851 36 : ALLOCATE (dcdr_env%matrix_m_alpha(3, nspins))
852 16 : DO i = 1, 3
853 34 : DO ispin = 1, nspins
854 18 : CALL cp_fm_create(dcdr_env%matrix_m_alpha(i, ispin), dcdr_env%likemos_fm_struct(1)%struct)
855 30 : CALL cp_fm_set_all(dcdr_env%matrix_m_alpha(i, ispin), 0.0_dp)
856 : END DO
857 : END DO
858 : END IF
859 :
860 : ! DBCSR matrices
861 24 : NULLIFY (dcdr_env%hamiltonian1)
862 24 : NULLIFY (dcdr_env%moments)
863 24 : NULLIFY (dcdr_env%matrix_difdip)
864 24 : NULLIFY (dcdr_env%matrix_core_charge_1)
865 24 : NULLIFY (dcdr_env%matrix_s1)
866 24 : NULLIFY (dcdr_env%matrix_t1)
867 24 : NULLIFY (dcdr_env%matrix_apply_op_constant)
868 24 : NULLIFY (dcdr_env%matrix_d_vhxc_dR)
869 24 : NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis)
870 24 : NULLIFY (dcdr_env%matrix_hc)
871 24 : NULLIFY (dcdr_env%matrix_ppnl_1)
872 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%perturbed_dm_correction, nspins)
873 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins)
874 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3)
875 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3)
876 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3)
877 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3)
878 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp2, 3)
879 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4)
880 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4)
881 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, nspins)
882 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3, nspins)
883 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6)
884 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3)
885 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3)
886 :
887 : ! temporary no_symmetry matrix:
888 96 : DO i = 1, 3
889 72 : CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp(i)%matrix)
890 : CALL dbcsr_create(dcdr_env%matrix_nosym_temp(i)%matrix, template=matrix_ks(1)%matrix, &
891 72 : matrix_type=dbcsr_type_no_symmetry)
892 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all)
893 72 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
894 :
895 72 : CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp2(i)%matrix)
896 : CALL dbcsr_create(dcdr_env%matrix_nosym_temp2(i)%matrix, template=matrix_ks(1)%matrix, &
897 72 : matrix_type=dbcsr_type_no_symmetry)
898 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp2(i)%matrix, sab_all)
899 96 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp2(i)%matrix, 0._dp)
900 : END DO
901 :
902 : ! moments carry the result of build_local_moment_matrix
903 96 : DO i = 1, 3
904 72 : CALL dbcsr_init_p(dcdr_env%moments(i)%matrix)
905 72 : CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments")
906 96 : CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp)
907 : END DO
908 24 : CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp])
909 :
910 96 : DO i = 1, 3
911 312 : DO j = 1, 3
912 216 : CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix)
913 216 : CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
914 288 : CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp)
915 : END DO
916 : END DO
917 :
918 52 : DO ispin = 1, nspins
919 28 : CALL dbcsr_init_p(dcdr_env%hamiltonian1(ispin)%matrix)
920 28 : CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction(ispin)%matrix)
921 28 : CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(ispin)%matrix)
922 28 : CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(ispin)%matrix, matrix_ks(1)%matrix)
923 52 : CALL dbcsr_copy(dcdr_env%perturbed_dm_correction(ispin)%matrix, matrix_ks(1)%matrix)
924 : END DO
925 :
926 : ! overlap/kinetic matrix: s(1) normal overlap matrix;
927 : ! s(2:4) derivatives wrt. nuclear coordinates
928 24 : CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix)
929 24 : CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix)
930 :
931 24 : CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix)
932 24 : CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix)
933 :
934 96 : DO i = 2, 4
935 72 : CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix)
936 72 : CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
937 :
938 72 : CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix)
939 96 : CALL dbcsr_copy(dcdr_env%matrix_t1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
940 : END DO
941 :
942 : ! j=1...3: derivative wrt nucleus A, 4...6: wrt nucleus B
943 52 : DO ispin = 1, nspins
944 220 : DO j = 1, 6
945 168 : CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix)
946 196 : CALL dbcsr_copy(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix, dcdr_env%matrix_s1(1)%matrix)
947 : END DO
948 : END DO
949 :
950 96 : DO i = 1, 3
951 72 : CALL dbcsr_init_p(dcdr_env%matrix_hc(i)%matrix)
952 : CALL dbcsr_create(dcdr_env%matrix_hc(i)%matrix, template=matrix_ks(1)%matrix, &
953 72 : matrix_type=dbcsr_type_symmetric)
954 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb)
955 96 : CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
956 : END DO
957 :
958 96 : DO i = 1, 3
959 72 : CALL dbcsr_init_p(dcdr_env%matrix_ppnl_1(i)%matrix)
960 : CALL dbcsr_create(dcdr_env%matrix_ppnl_1(i)%matrix, template=matrix_ks(1)%matrix, &
961 72 : matrix_type=dbcsr_type_symmetric)
962 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb)
963 96 : CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
964 : END DO
965 :
966 96 : DO i = 1, 3
967 156 : DO ispin = 1, nspins
968 84 : CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix)
969 156 : CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, dcdr_env%matrix_s1(1)%matrix)
970 : END DO
971 :
972 72 : CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix)
973 72 : CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix)
974 96 : CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
975 : END DO
976 :
977 : ! CHC
978 52 : DO ispin = 1, nspins
979 28 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
980 28 : CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
981 :
982 28 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, nmo)
983 : ! chc = mo * matrix_ks * mo
984 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, &
985 : 1.0_dp, mo_coeff, buf, &
986 28 : 0.0_dp, dcdr_env%chc(ispin))
987 :
988 80 : CALL cp_fm_release(buf)
989 : END DO
990 :
991 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
992 24 : "PRINT%PROGRAM_RUN_INFO")
993 :
994 24 : CALL timestop(handle)
995 72 : END SUBROUTINE dcdr_env_init
996 :
997 : ! **************************************************************************************************
998 : !> \brief Deallocate the dcdr environment
999 : !> \param qs_env ...
1000 : !> \param dcdr_env ...
1001 : ! **************************************************************************************************
1002 24 : SUBROUTINE dcdr_env_cleanup(qs_env, dcdr_env)
1003 :
1004 : TYPE(qs_environment_type), POINTER :: qs_env
1005 : TYPE(dcdr_env_type) :: dcdr_env
1006 :
1007 : INTEGER :: ispin
1008 : TYPE(cp_logger_type), POINTER :: logger
1009 : TYPE(section_vals_type), POINTER :: dcdr_section
1010 :
1011 : ! Destroy the logger
1012 24 : logger => cp_get_default_logger()
1013 24 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
1014 24 : CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT")
1015 :
1016 24 : DEALLOCATE (dcdr_env%list_of_atoms)
1017 :
1018 24 : CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct)
1019 24 : CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct)
1020 52 : DO ispin = 1, dcdr_env%nspins
1021 28 : CALL cp_fm_struct_release(dcdr_env%momo_fm_struct(ispin)%struct)
1022 52 : CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct(ispin)%struct)
1023 : END DO
1024 24 : DEALLOCATE (dcdr_env%momo_fm_struct)
1025 24 : DEALLOCATE (dcdr_env%likemos_fm_struct)
1026 :
1027 24 : DEALLOCATE (dcdr_env%deltar)
1028 24 : DEALLOCATE (dcdr_env%delta_basis_function)
1029 :
1030 24 : IF (dcdr_env%localized_psi0) THEN
1031 : ! DEALLOCATE (dcdr_env%psi0_order)
1032 4 : DEALLOCATE (dcdr_env%centers_set(1)%array)
1033 4 : DEALLOCATE (dcdr_env%center_list(1)%array)
1034 4 : DEALLOCATE (dcdr_env%centers_set)
1035 4 : DEALLOCATE (dcdr_env%center_list)
1036 4 : DEALLOCATE (dcdr_env%apt_subset)
1037 : END IF
1038 :
1039 24 : DEALLOCATE (dcdr_env%apt_el_dcdr)
1040 24 : DEALLOCATE (dcdr_env%apt_nuc_dcdr)
1041 24 : DEALLOCATE (dcdr_env%apt_total_dcdr)
1042 24 : IF (dcdr_env%localized_psi0) THEN
1043 4 : DEALLOCATE (dcdr_env%apt_el_dcdr_per_center)
1044 4 : DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset)
1045 : END IF
1046 :
1047 : ! Full matrices
1048 24 : CALL cp_fm_release(dcdr_env%dCR)
1049 24 : CALL cp_fm_release(dcdr_env%dCR_prime)
1050 24 : CALL cp_fm_release(dcdr_env%mo_coeff)
1051 24 : CALL cp_fm_release(dcdr_env%chc)
1052 24 : CALL cp_fm_release(dcdr_env%op_dR)
1053 :
1054 24 : IF (dcdr_env%z_matrix_method) THEN
1055 4 : CALL cp_fm_release(dcdr_env%matrix_m_alpha)
1056 : END IF
1057 :
1058 : ! DBCSR matrices
1059 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
1060 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)
1061 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%moments)
1062 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip)
1063 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1)
1064 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp)
1065 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp2)
1066 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s)
1067 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t)
1068 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1)
1069 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1)
1070 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant)
1071 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR)
1072 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis)
1073 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc)
1074 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1)
1075 :
1076 24 : END SUBROUTINE dcdr_env_cleanup
1077 :
1078 : END MODULE qs_dcdr_utils
|
