(git:374b731)
Loading...
Searching...
No Matches
qs_dcdr_utils.F
Go to the documentation of this file.
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
13MODULE 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
19 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
20 dbcsr_allocate_matrix_set,&
21 dbcsr_deallocate_matrix_set
22 USE cp_files, ONLY: close_file,&
23 open_file
24 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 cp_fm_struct_release
26 USE cp_fm_types, ONLY: cp_fm_create,&
27 cp_fm_get_info,&
28 cp_fm_get_submatrix,&
29 cp_fm_release,&
30 cp_fm_set_all,&
31 cp_fm_set_submatrix,&
32 cp_fm_to_fm,&
33 cp_fm_type
34 USE cp_log_handling, ONLY: cp_get_default_logger,&
35 cp_logger_get_default_io_unit,&
36 cp_logger_type,&
37 cp_to_string
38 USE cp_output_handling, ONLY: cp_p_file,&
39 cp_print_key_finished_output,&
40 cp_print_key_generate_filename,&
41 cp_print_key_should_output,&
42 cp_print_key_unit_nr
43 USE cp_result_methods, ONLY: get_results
44 USE cp_result_types, ONLY: cp_result_type
45 USE dbcsr_api, ONLY: dbcsr_copy,&
46 dbcsr_create,&
47 dbcsr_init_p,&
48 dbcsr_p_type,&
49 dbcsr_set,&
50 dbcsr_type,&
51 dbcsr_type_no_symmetry,&
52 dbcsr_type_symmetric
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
100CONTAINS
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 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, work
115 INTEGER, INTENT(IN) :: nmo, icenter
116 TYPE(cp_fm_type), INTENT(IN) :: res
117
118 CHARACTER(LEN=*), PARAMETER :: routinen = 'multiply_localization'
119
120 INTEGER :: handle
121
122 CALL timeset(routinen, handle)
123
124 ! Multiply by the MO coefficients
125 CALL cp_dbcsr_sm_fm_multiply(ao_matrix, mo_coeff, work, ncol=nmo)
126
127 ! Only keep the icenter-th column
128 CALL cp_fm_to_fm(work, res, 1, icenter, icenter)
129
130 ! Reset the matrices
131 CALL cp_fm_set_all(work, 0.0_dp)
132
133 CALL timestop(handle)
134 END SUBROUTINE multiply_localization
135
136! **************************************************************************************************
137!> \brief Copied from linres_read_restart
138!> \param qs_env ...
139!> \param linres_section ...
140!> \param vec ...
141!> \param lambda ...
142!> \param beta ...
143!> \param tag ...
144!> \note Adapted from linres_read_restart (ED)
145!> Would be nice not to crash but to start from zero if the present file doesn't match.
146! **************************************************************************************************
147 SUBROUTINE dcdr_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
148 TYPE(qs_environment_type), POINTER :: qs_env
149 TYPE(section_vals_type), POINTER :: linres_section
150 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
151 INTEGER, INTENT(IN) :: lambda, beta
152 CHARACTER(LEN=*) :: tag
153
154 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_read_restart'
155
156 CHARACTER(LEN=default_path_length) :: filename
157 CHARACTER(LEN=default_string_length) :: my_middle
158 INTEGER :: beta_tmp, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, lambda_tmp, &
159 max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
160 LOGICAL :: file_exists
161 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
162 TYPE(cp_fm_type), POINTER :: mo_coeff
163 TYPE(cp_logger_type), POINTER :: logger
164 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
165 TYPE(mp_para_env_type), POINTER :: para_env
166 TYPE(section_vals_type), POINTER :: print_key
167
168 file_exists = .false.
169
170 CALL timeset(routinen, handle)
171
172 NULLIFY (mos, para_env, logger, print_key, vecbuffer)
173 logger => cp_get_default_logger()
174
175 iounit = cp_print_key_unit_nr(logger, linres_section, &
176 "PRINT%PROGRAM_RUN_INFO", extension=".Log")
177
178 CALL get_qs_env(qs_env=qs_env, &
179 para_env=para_env, &
180 mos=mos)
181
182 nspins = SIZE(mos)
183
184 rst_unit = -1
185 IF (para_env%is_source()) THEN
186 CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
187 n_rep_val=n_rep_val)
188
189 CALL xstring(tag, ia, ie)
190 my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
191 //trim("-")//trim(adjustl(cp_to_string(lambda)))
192
193 IF (n_rep_val > 0) THEN
194 CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
195 CALL xstring(filename, ia, ie)
196 filename = filename(ia:ie)//trim(my_middle)//".lr"
197 ELSE
198 ! try to read from the filename that is generated automatically from the printkey
199 print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
200 filename = cp_print_key_generate_filename(logger, print_key, &
201 extension=".lr", middle_name=trim(my_middle), my_local=.false.)
202 END IF
203 INQUIRE (file=filename, exist=file_exists)
204 !
205 ! open file
206 IF (file_exists) THEN
207 CALL open_file(file_name=trim(filename), &
208 file_action="READ", &
209 file_form="UNFORMATTED", &
210 file_position="REWIND", &
211 file_status="OLD", &
212 unit_number=rst_unit)
213
214 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
215 "LINRES| Reading response wavefunctions from the restart file <"//trim(adjustl(filename))//">"
216 ELSE
217 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
218 "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
219 END IF
220 END IF
221
222 CALL para_env%bcast(file_exists)
223
224 IF (file_exists) THEN
225
226 CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
227 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
228
229 ALLOCATE (vecbuffer(nao, max_block))
230 !
231 ! read headers
232 IF (rst_unit > 0) READ (rst_unit, iostat=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
233 CALL para_env%bcast(iostat)
234 CALL para_env%bcast(beta_tmp)
235 CALL para_env%bcast(lambda_tmp)
236 CALL para_env%bcast(nspins_tmp)
237 CALL para_env%bcast(nao_tmp)
238
239 ! check that the number nao, nmo and nspins are
240 ! the same as in the current mos
241 IF (nspins_tmp .NE. nspins) THEN
242 cpabort("nspins not consistent")
243 END IF
244 IF (nao_tmp .NE. nao) cpabort("nao not consistent")
245 ! check that it's the right file
246 ! the same as in the current mos
247 IF (lambda_tmp .NE. lambda) cpabort("lambda not consistent")
248 IF (beta_tmp .NE. beta) cpabort("beta not consistent")
249 !
250 DO ispin = 1, nspins
251 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
252 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
253 !
254 IF (rst_unit > 0) READ (rst_unit) nmo_tmp
255 CALL para_env%bcast(nmo_tmp)
256 IF (nmo_tmp .NE. nmo) cpabort("nmo not consistent")
257 !
258 ! read the response
259 DO i = 1, nmo, max(max_block, 1)
260 i_block = min(max_block, nmo - i + 1)
261 DO j = 1, i_block
262 IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
263 END DO
264 CALL para_env%bcast(vecbuffer)
265 CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
266 END DO
267 END DO
268
269 IF (iostat /= 0) THEN
270 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
271 "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
272 END IF
273
274 DEALLOCATE (vecbuffer)
275
276 END IF
277
278 IF (para_env%is_source()) THEN
279 IF (file_exists) CALL close_file(unit_number=rst_unit)
280 END IF
281
282 CALL timestop(handle)
283
284 END SUBROUTINE dcdr_read_restart
285
286! **************************************************************************************************
287!> \brief Copied from linres_write_restart
288!> \param qs_env ...
289!> \param linres_section ...
290!> \param vec ...
291!> \param lambda ...
292!> \param beta ...
293!> \param tag ...
294!> \note Adapted from linres_read_restart (ED)
295!> Would be nice not to crash but to start from zero if the present file doesn't match.
296! **************************************************************************************************
297 SUBROUTINE dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
298 TYPE(qs_environment_type), POINTER :: qs_env
299 TYPE(section_vals_type), POINTER :: linres_section
300 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
301 INTEGER, INTENT(IN) :: lambda, beta
302 CHARACTER(LEN=*) :: tag
303
304 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_write_restart'
305
306 CHARACTER(LEN=default_path_length) :: filename
307 CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status
308 INTEGER :: handle, i, i_block, ia, ie, iounit, &
309 ispin, j, max_block, nao, nmo, nspins, &
310 rst_unit
311 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
312 TYPE(cp_fm_type), POINTER :: mo_coeff
313 TYPE(cp_logger_type), POINTER :: logger
314 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
315 TYPE(mp_para_env_type), POINTER :: para_env
316 TYPE(section_vals_type), POINTER :: print_key
317
318 NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
319
320 CALL timeset(routinen, handle)
321
322 logger => cp_get_default_logger()
323
324 IF (btest(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
325 used_print_key=print_key), &
326 cp_p_file)) THEN
327
328 iounit = cp_print_key_unit_nr(logger, linres_section, &
329 "PRINT%PROGRAM_RUN_INFO", extension=".Log")
330
331 CALL get_qs_env(qs_env=qs_env, &
332 mos=mos, &
333 para_env=para_env)
334
335 nspins = SIZE(mos)
336
337 my_status = "REPLACE"
338 my_pos = "REWIND"
339 CALL xstring(tag, ia, ie)
340 my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
341 //trim("-")//trim(adjustl(cp_to_string(lambda)))
342 rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
343 extension=".lr", middle_name=trim(my_middle), file_status=trim(my_status), &
344 file_position=trim(my_pos), file_action="WRITE", file_form="UNFORMATTED")
345
346 filename = cp_print_key_generate_filename(logger, print_key, &
347 extension=".lr", middle_name=trim(my_middle), my_local=.false.)
348
349 IF (iounit > 0) THEN
350 WRITE (unit=iounit, fmt="(T2,A)") &
351 "LINRES| Writing response functions to the restart file <"//trim(adjustl(filename))//">"
352 END IF
353
354 !
355 ! write data to file
356 ! use the scalapack block size as a default for buffering columns
357 CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
358 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
359 ALLOCATE (vecbuffer(nao, max_block))
360
361 IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
362
363 DO ispin = 1, nspins
364 CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
365
366 IF (rst_unit > 0) WRITE (rst_unit) nmo
367
368 DO i = 1, nmo, max(max_block, 1)
369 i_block = min(max_block, nmo - i + 1)
370 CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
371 ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
372 ! to old ones, and in cases where max_block is different between runs, as might happen during
373 ! restarts with a different number of CPUs
374 DO j = 1, i_block
375 IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
376 END DO
377 END DO
378 END DO
379
380 DEALLOCATE (vecbuffer)
381
382 CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
383 "PRINT%RESTART")
384 END IF
385
386 CALL timestop(handle)
387
388 END SUBROUTINE dcdr_write_restart
389
390! **************************************************************************************************
391!> \brief Print the APT and sum rules
392!> \param dcdr_env ...
393!> \param qs_env ...
394!> \author Edward Ditler
395! **************************************************************************************************
396 SUBROUTINE dcdr_print(dcdr_env, qs_env)
397 TYPE(dcdr_env_type) :: dcdr_env
398 TYPE(qs_environment_type), POINTER :: qs_env
399
400 CHARACTER(LEN=default_string_length) :: description
401 INTEGER :: alpha, beta, delta, gamma, i, k, l, &
402 lambda, natom, nsubset, output_unit
403 REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
404 REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center_dcdr, apt_subset_dcdr
405 REAL(kind=dp), DIMENSION(3, 3) :: sum_rule_0, sum_rule_1, sum_rule_2
406 TYPE(cp_logger_type), POINTER :: logger
407 TYPE(cp_result_type), POINTER :: results
408 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
409 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
410 TYPE(section_vals_type), POINTER :: dcdr_section
411
412 NULLIFY (logger)
413
414 logger => cp_get_default_logger()
415 output_unit = cp_logger_get_default_io_unit(logger)
416
417 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
418
419 NULLIFY (particle_set)
420 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
421 natom = SIZE(particle_set)
422 nsubset = SIZE(molecule_set)
423
424 apt_el_dcdr => dcdr_env%apt_el_dcdr
425 apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr
426 apt_total_dcdr => dcdr_env%apt_total_dcdr
427 apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset
428 apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center
429
430 IF (dcdr_env%localized_psi0) THEN
431 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the final apt matrix per atom per subset'
432 DO k = 1, natom
433 DO l = 1, nsubset
434 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, *) 'APT | Subset', l
435 DO i = 1, 3
436 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,F15.6,F15.6,F15.6)") &
437 'APT | apt_subset ', i, apt_subset_dcdr(i, :, k, l)
438 END DO
439 END DO
440 END DO
441 END IF
442
443 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") &
444 'APT | Write the final apt matrix per atom (Position perturbation)'
445 DO l = 1, natom
446 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") &
447 'APT | Atom', l, ' - GAPT ', &
448 (apt_total_dcdr(1, 1, l) &
449 + apt_total_dcdr(2, 2, l) &
450 + apt_total_dcdr(3, 3, l))/3._dp
451 DO i = 1, 3
452 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
453 END DO
454 END DO
455
456 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix'
457 DO i = 1, 3
458 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
459 "(A,F15.6,F15.6,F15.6)") "APT | ", sum(apt_total_dcdr(i, :, :), dim=2)
460 END DO
461 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | End Write the final apt matrix'
462
463 ! Get the dipole
464 CALL get_qs_env(qs_env, results=results)
465 description = "[DIPOLE]"
466 CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3))
467
468 ! Sum rules [for all alpha, beta]
469 sum_rule_0 = 0._dp
470 sum_rule_1 = 0._dp
471 sum_rule_2 = 0._dp
472
473 DO alpha = 1, 3
474 DO beta = 1, 3
475 ! 0: sum_lambda apt(alpha, beta, lambda)
476 DO lambda = 1, natom
477 sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
478 + apt_total_dcdr(alpha, beta, lambda)
479 END DO
480
481 ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
482 DO gamma = 1, 3
483 sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
484 + levi_civita(alpha, beta, gamma)*dcdr_env%dipole_pos(gamma)
485 END DO
486
487 ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
488 DO lambda = 1, natom
489 DO gamma = 1, 3
490 DO delta = 1, 3
491 sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
492 + levi_civita(beta, gamma, delta) &
493 *particle_set(lambda)%r(gamma) &
494 *apt_total_dcdr(delta, alpha, lambda)
495 END DO
496 END DO
497 END DO
498
499 END DO ! beta
500 END DO ! alpha
501
502 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules"
503 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") &
504 "APT |", "Total APT", "Dipole", "R * APT"
505 DO alpha = 1, 3
506 DO beta = 1, 3
507 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
508 "(A,I3,I3,F15.6,F15.6,F15.6)") &
509 "APT | ", &
510 alpha, beta, &
511 sum_rule_0(alpha, beta), &
512 sum_rule_1(alpha, beta), &
513 sum_rule_2(alpha, beta)
514 END DO
515 END DO
516
517 END SUBROUTINE dcdr_print
518
519! **************************************************************************************************
520!> \brief ...
521!> \param r ...
522!> \param cell ...
523!> \param r_shifted ...
524! **************************************************************************************************
525 SUBROUTINE shift_wannier_into_cell(r, cell, r_shifted)
526 REAL(dp), DIMENSION(3), INTENT(in) :: r
527 TYPE(cell_type), INTENT(in), POINTER :: cell
528 REAL(dp), DIMENSION(3), INTENT(out) :: r_shifted
529
530 INTEGER :: i
531 REAL(kind=dp), DIMENSION(3) :: abc
532
533 ! Only orthorombic cell for now
534 CALL get_cell(cell, abc=abc)
535
536 DO i = 1, 3
537 IF (r(i) < 0._dp) THEN
538 r_shifted(i) = r(i) + abc(i)
539 ELSE IF (r(i) > abc(i)) THEN
540 r_shifted(i) = r(i) - abc(i)
541 ELSE
542 r_shifted(i) = r(i)
543 END IF
544 END DO
545 END SUBROUTINE shift_wannier_into_cell
546
547! **************************************************************************************************
548!> \brief ...
549!> \param dcdr_env ...
550!> \param qs_env ...
551! **************************************************************************************************
552 SUBROUTINE get_loc_setting(dcdr_env, qs_env)
553 TYPE(dcdr_env_type) :: dcdr_env
554 TYPE(qs_environment_type), POINTER :: qs_env
555
556 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_loc_setting'
557
558 INTEGER :: handle, is, ispin, istate, max_states, &
559 nmo, nmoloc, nstate, nstate_list(2)
560 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
561 REAL(dp), DIMENSION(:, :), POINTER :: center_array
562 TYPE(linres_control_type), POINTER :: linres_control
563 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
564 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
565 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
566 TYPE(section_vals_type), POINTER :: dcdr_section
567
568 CALL timeset(routinen, handle)
569
570 CALL get_qs_env(qs_env=qs_env, &
571 linres_control=linres_control, &
572 mos=mos)
573
574 ! Some checks
575 max_states = 0
576 CALL get_mo_set(mo_set=mos(1), nmo=nmo)
577 max_states = max(max_states, nmo)
578
579 ! check that the number of localized states is equal to the number of states
580 nmoloc = SIZE(linres_control%qs_loc_env%localized_wfn_control%centers_set(1)%array, 2)
581 IF (nmoloc .NE. nmo) THEN
582 cpabort("The number of localized functions is not equal to the number of states.")
583 END IF
584
585 ! which center for the orbitals shall we use
586 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
587 CALL section_vals_val_get(dcdr_section, "ORBITAL_CENTER", i_val=dcdr_env%orb_center)
588 SELECT CASE (dcdr_env%orb_center)
589 CASE (current_orb_center_wannier)
590 dcdr_env%orb_center_name = "WANNIER"
591 CASE DEFAULT
592 cpabort(" ")
593 END SELECT
594
595 qs_loc_env => linres_control%qs_loc_env
596 CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
597
598 ALLOCATE (dcdr_env%centers_set(dcdr_env%nspins))
599 ALLOCATE (dcdr_env%center_list(dcdr_env%nspins))
600 ALLOCATE (state_list(max_states, dcdr_env%nspins))
601 state_list(:, :) = huge(0)
602 nstate_list(:) = huge(0)
603
604 ! Build the state_list
605 DO ispin = 1, dcdr_env%nspins
606 center_array => localized_wfn_control%centers_set(ispin)%array
607 nstate = 0
608 DO istate = 1, SIZE(center_array, 2)
609 nstate = nstate + 1
610 state_list(nstate, ispin) = istate
611 END DO
612 nstate_list(ispin) = nstate
613
614 ! clustering the states
615 nstate = nstate_list(ispin)
616 dcdr_env%nstates(ispin) = nstate
617
618 ALLOCATE (dcdr_env%center_list(ispin)%array(2, nstate + 1))
619 ALLOCATE (dcdr_env%centers_set(ispin)%array(3, nstate))
620 dcdr_env%center_list(ispin)%array(:, :) = huge(0)
621 dcdr_env%centers_set(ispin)%array(:, :) = huge(0.0_dp)
622
623 center_array => localized_wfn_control%centers_set(ispin)%array
624
625 ! point to the psi0 centers
626 SELECT CASE (dcdr_env%orb_center)
627 CASE (current_orb_center_wannier)
628 ! use the wannier center as -center-
629 dcdr_env%nbr_center(ispin) = nstate
630 DO is = 1, nstate
631 istate = state_list(is, 1)
632 dcdr_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
633 dcdr_env%center_list(ispin)%array(1, is) = is
634 dcdr_env%center_list(ispin)%array(2, is) = istate
635 END DO
636 dcdr_env%center_list(ispin)%array(1, nstate + 1) = nstate + 1
637
638 CASE DEFAULT
639 cpabort("Unknown orbital center...")
640 END SELECT
641 END DO
642
643 CALL timestop(handle)
644 END SUBROUTINE get_loc_setting
645
646! **************************************************************************************************
647!> \brief Initialize the dcdr environment
648!> \param dcdr_env ...
649!> \param qs_env ...
650! **************************************************************************************************
651 SUBROUTINE dcdr_env_init(dcdr_env, qs_env)
652 TYPE(dcdr_env_type) :: dcdr_env
653 TYPE(qs_environment_type), POINTER :: qs_env
654
655 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_env_init'
656
657 INTEGER :: handle, homo, i, isize, ispin, j, jg, &
658 n_rep, nao, natom, nmo, nspins, &
659 nsubset, output_unit, reference, &
660 unit_number
661 INTEGER, DIMENSION(:), POINTER :: tmplist
662 LOGICAL :: explicit
663 REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
664 TYPE(cp_fm_type) :: buf
665 TYPE(cp_fm_type), POINTER :: mo_coeff
666 TYPE(cp_logger_type), POINTER :: logger
667 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
668 TYPE(dft_control_type), POINTER :: dft_control
669 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
670 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
671 TYPE(mp_para_env_type), POINTER :: para_env
672 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
673 POINTER :: sab_all, sab_orb
674 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
675 TYPE(qs_ks_env_type), POINTER :: ks_env
676 TYPE(section_vals_type), POINTER :: dcdr_section, loc_section, lr_section
677
678 CALL timeset(routinen, handle)
679
680 ! Set up the logger
681 NULLIFY (logger, loc_section, dcdr_section, lr_section)
682 logger => cp_get_default_logger()
683 loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE")
684 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
685 dcdr_env%output_unit = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", &
686 extension=".data", middle_name="dcdr", log_filename=.false., &
687 file_position="REWIND", file_status="REPLACE")
688
689 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
690 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
691 extension=".linresLog")
692 unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
693
694 IF (output_unit > 0) THEN
695 WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***"
696 END IF
697
698 NULLIFY (ks_env, dft_control, sab_orb, sab_all, particle_set, molecule_set, matrix_s, matrix_ks, mos, para_env)
699 CALL get_qs_env(qs_env=qs_env, &
700 ks_env=ks_env, &
701 dft_control=dft_control, &
702 sab_orb=sab_orb, &
703 sab_all=sab_all, &
704 particle_set=particle_set, &
705 molecule_set=molecule_set, &
706 matrix_s=matrix_s, &
707 matrix_ks=matrix_ks, &
708 mos=mos, &
709 para_env=para_env)
710
711 natom = SIZE(particle_set)
712 nsubset = SIZE(molecule_set)
713 nspins = dft_control%nspins
714 dcdr_env%nspins = dft_control%nspins
715
716 NULLIFY (dcdr_env%matrix_s)
717 CALL build_overlap_matrix(ks_env, matrix_s=dcdr_env%matrix_s, &
718 matrix_name="OVERLAP MATRIX", &
719 nderivative=1, &
720 basis_type_a="ORB", &
721 basis_type_b="ORB", &
722 sab_nl=sab_orb)
723
724 NULLIFY (dcdr_env%matrix_t)
725 CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, &
726 matrix_name="KINETIC ENERGY MATRIX", &
727 basis_type="ORB", &
728 sab_nl=sab_orb, nderivative=1, &
729 eps_filter=dft_control%qs_control%eps_filter_matrix)
730
731 ! Get inputs
732 CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
733 CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
734 CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
735 dcdr_env%ref_point = 0._dp
736
737 ! List of atoms
738 NULLIFY (tmplist)
739 isize = 0
740 CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep)
741 IF (n_rep == 0) THEN
742 ALLOCATE (dcdr_env%list_of_atoms(natom))
743 DO jg = 1, natom
744 dcdr_env%list_of_atoms(jg) = jg
745 END DO
746 ELSE
747 DO jg = 1, n_rep
748 ALLOCATE (dcdr_env%list_of_atoms(isize))
749 CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist)
750 CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist))
751 dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
752 isize = SIZE(dcdr_env%list_of_atoms)
753 END DO
754 END IF
755
756 ! Reference point
757 IF (dcdr_env%localized_psi0) THEN
758 ! Get the Wannier localized wave functions and centers
759 CALL get_loc_setting(dcdr_env, qs_env)
760 ELSE
761 ! Get the reference point from the input
762 CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference)
763 CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit)
764 IF (explicit) THEN
765 CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point)
766 ELSE
767 IF (reference == use_mom_ref_user) &
768 cpabort("User-defined reference point should be given explicitly")
769 END IF
770
771 CALL get_reference_point(rpoint=dcdr_env%ref_point, qs_env=qs_env, &
772 reference=reference, &
773 ref_point=ref_point)
774 END IF
775
776 ! Helper matrix structs
777 NULLIFY (dcdr_env%aoao_fm_struct, &
778 dcdr_env%momo_fm_struct, &
779 dcdr_env%likemos_fm_struct, &
780 dcdr_env%homohomo_fm_struct)
781 CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, &
782 nao=nao, nmo=nmo, homo=homo)
783 CALL cp_fm_struct_create(dcdr_env%aoao_fm_struct, nrow_global=nao, &
784 ncol_global=nao, para_env=para_env, &
785 context=mo_coeff%matrix_struct%context)
786 CALL cp_fm_struct_create(dcdr_env%homohomo_fm_struct, nrow_global=homo, &
787 ncol_global=homo, para_env=para_env, &
788 context=mo_coeff%matrix_struct%context)
789 dcdr_env%nao = nao
790
791 ALLOCATE (dcdr_env%nmo(nspins))
792 ALLOCATE (dcdr_env%momo_fm_struct(nspins))
793 ALLOCATE (dcdr_env%likemos_fm_struct(nspins))
794 DO ispin = 1, nspins
795 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
796 nao=nao, nmo=nmo, homo=homo)
797 CALL cp_fm_struct_create(dcdr_env%momo_fm_struct(ispin)%struct, nrow_global=nmo, &
798 ncol_global=nmo, para_env=para_env, &
799 context=mo_coeff%matrix_struct%context)
800 CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct(ispin)%struct, &
801 template_fmstruct=mo_coeff%matrix_struct)
802 dcdr_env%nmo(ispin) = nmo
803 END DO
804
805 ! Fields of reals
806 ALLOCATE (dcdr_env%deltaR(3, natom))
807 ALLOCATE (dcdr_env%delta_basis_function(3, natom))
808 ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom))
809 ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom))
810 ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom))
811
812 dcdr_env%apt_el_dcdr = 0._dp
813 dcdr_env%apt_nuc_dcdr = 0._dp
814 dcdr_env%apt_total_dcdr = 0._dp
815
816 dcdr_env%deltaR = 0.0_dp
817 dcdr_env%delta_basis_function = 0._dp
818
819 ! Localization
820 IF (dcdr_env%localized_psi0) THEN
821 ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1)))
822 ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset))
823 ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset))
824 dcdr_env%apt_el_dcdr_per_center = 0._dp
825 dcdr_env%apt_el_dcdr_per_subset = 0._dp
826 dcdr_env%apt_subset = 0.0_dp
827 END IF
828
829 ! Full matrices
830 ALLOCATE (dcdr_env%mo_coeff(nspins))
831 ALLOCATE (dcdr_env%dCR(nspins))
832 ALLOCATE (dcdr_env%dCR_prime(nspins))
833 ALLOCATE (dcdr_env%chc(nspins))
834 ALLOCATE (dcdr_env%op_dR(nspins))
835
836 DO ispin = 1, nspins
837 CALL cp_fm_create(dcdr_env%dCR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
838 CALL cp_fm_create(dcdr_env%dCR_prime(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
839 CALL cp_fm_create(dcdr_env%mo_coeff(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
840 CALL cp_fm_create(dcdr_env%chc(ispin), dcdr_env%momo_fm_struct(ispin)%struct)
841 CALL cp_fm_create(dcdr_env%op_dR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
842
843 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
844 CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
845 END DO
846
847 ! DBCSR matrices
848 NULLIFY (dcdr_env%hamiltonian1)
849 NULLIFY (dcdr_env%moments)
850 NULLIFY (dcdr_env%matrix_difdip)
851 NULLIFY (dcdr_env%matrix_core_charge_1)
852 NULLIFY (dcdr_env%matrix_s1)
853 NULLIFY (dcdr_env%matrix_t1)
854 NULLIFY (dcdr_env%matrix_apply_op_constant)
855 NULLIFY (dcdr_env%matrix_d_vhxc_dR)
856 NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis)
857 NULLIFY (dcdr_env%matrix_hc)
858 NULLIFY (dcdr_env%matrix_ppnl_1)
859 CALL dbcsr_allocate_matrix_set(dcdr_env%perturbed_dm_correction, nspins)
860 CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins)
861 CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3)
862 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3)
863 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3)
864 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3)
865 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp2, 3)
866 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4)
867 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4)
868 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, nspins)
869 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3, nspins)
870 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6)
871 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3)
872 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3)
873
874 ! temporary no_symmetry matrix:
875 DO i = 1, 3
876 CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp(i)%matrix)
877 CALL dbcsr_create(dcdr_env%matrix_nosym_temp(i)%matrix, template=matrix_ks(1)%matrix, &
878 matrix_type=dbcsr_type_no_symmetry)
879 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all)
880 CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
881
882 CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp2(i)%matrix)
883 CALL dbcsr_create(dcdr_env%matrix_nosym_temp2(i)%matrix, template=matrix_ks(1)%matrix, &
884 matrix_type=dbcsr_type_no_symmetry)
885 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp2(i)%matrix, sab_all)
886 CALL dbcsr_set(dcdr_env%matrix_nosym_temp2(i)%matrix, 0._dp)
887 END DO
888
889 ! moments carry the result of build_local_moment_matrix
890 DO i = 1, 3
891 CALL dbcsr_init_p(dcdr_env%moments(i)%matrix)
892 CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments")
893 CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp)
894 END DO
895 CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp])
896
897 DO i = 1, 3
898 DO j = 1, 3
899 CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix)
900 CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
901 CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp)
902 END DO
903 END DO
904
905 DO ispin = 1, nspins
906 CALL dbcsr_init_p(dcdr_env%hamiltonian1(ispin)%matrix)
907 CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction(ispin)%matrix)
908 CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(ispin)%matrix)
909 CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(ispin)%matrix, matrix_ks(1)%matrix)
910 CALL dbcsr_copy(dcdr_env%perturbed_dm_correction(ispin)%matrix, matrix_ks(1)%matrix)
911 END DO
912
913 ! overlap/kinetic matrix: s(1) normal overlap matrix;
914 ! s(2:4) derivatives wrt. nuclear coordinates
915 CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix)
916 CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix)
917
918 CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix)
919 CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix)
920
921 DO i = 2, 4
922 CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix)
923 CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
924
925 CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix)
926 CALL dbcsr_copy(dcdr_env%matrix_t1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
927 END DO
928
929 ! j=1...3: derivative wrt nucleus A, 4...6: wrt nucleus B
930 DO ispin = 1, nspins
931 DO j = 1, 6
932 CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix)
933 CALL dbcsr_copy(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix, dcdr_env%matrix_s1(1)%matrix)
934 END DO
935 END DO
936
937 DO i = 1, 3
938 CALL dbcsr_init_p(dcdr_env%matrix_hc(i)%matrix)
939 CALL dbcsr_create(dcdr_env%matrix_hc(i)%matrix, template=matrix_ks(1)%matrix, &
940 matrix_type=dbcsr_type_symmetric)
941 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb)
942 CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
943 END DO
944
945 DO i = 1, 3
946 CALL dbcsr_init_p(dcdr_env%matrix_ppnl_1(i)%matrix)
947 CALL dbcsr_create(dcdr_env%matrix_ppnl_1(i)%matrix, template=matrix_ks(1)%matrix, &
948 matrix_type=dbcsr_type_symmetric)
949 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb)
950 CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
951 END DO
952
953 DO i = 1, 3
954 DO ispin = 1, nspins
955 CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix)
956 CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, dcdr_env%matrix_s1(1)%matrix)
957 END DO
958
959 CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix)
960 CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix)
961 CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
962 END DO
963
964 ! CHC
965 DO ispin = 1, nspins
966 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
967 CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
968
969 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, nmo)
970 ! chc = mo * matrix_ks * mo
971 CALL parallel_gemm('T', 'N', nmo, nmo, nao, &
972 1.0_dp, mo_coeff, buf, &
973 0.0_dp, dcdr_env%chc(ispin))
974
975 CALL cp_fm_release(buf)
976 END DO
977
978 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
979 "PRINT%PROGRAM_RUN_INFO")
980
981 CALL timestop(handle)
982 END SUBROUTINE dcdr_env_init
983
984! **************************************************************************************************
985!> \brief Deallocate the dcdr environment
986!> \param qs_env ...
987!> \param dcdr_env ...
988! **************************************************************************************************
989 SUBROUTINE dcdr_env_cleanup(qs_env, dcdr_env)
990
991 TYPE(qs_environment_type), POINTER :: qs_env
992 TYPE(dcdr_env_type) :: dcdr_env
993
994 INTEGER :: ispin
995 TYPE(cp_logger_type), POINTER :: logger
996 TYPE(section_vals_type), POINTER :: dcdr_section
997
998 ! Destroy the logger
999 logger => cp_get_default_logger()
1000 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
1001 CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT")
1002
1003 DEALLOCATE (dcdr_env%list_of_atoms)
1004
1005 CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct)
1006 CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct)
1007 DO ispin = 1, dcdr_env%nspins
1008 CALL cp_fm_struct_release(dcdr_env%momo_fm_struct(ispin)%struct)
1009 CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct(ispin)%struct)
1010 END DO
1011 DEALLOCATE (dcdr_env%momo_fm_struct)
1012 DEALLOCATE (dcdr_env%likemos_fm_struct)
1013
1014 DEALLOCATE (dcdr_env%deltar)
1015 DEALLOCATE (dcdr_env%delta_basis_function)
1016
1017 IF (dcdr_env%localized_psi0) THEN
1018 ! DEALLOCATE (dcdr_env%psi0_order)
1019 DEALLOCATE (dcdr_env%centers_set(1)%array)
1020 DEALLOCATE (dcdr_env%center_list(1)%array)
1021 DEALLOCATE (dcdr_env%centers_set)
1022 DEALLOCATE (dcdr_env%center_list)
1023 DEALLOCATE (dcdr_env%apt_subset)
1024 END IF
1025
1026 DEALLOCATE (dcdr_env%apt_el_dcdr)
1027 DEALLOCATE (dcdr_env%apt_nuc_dcdr)
1028 DEALLOCATE (dcdr_env%apt_total_dcdr)
1029 IF (dcdr_env%localized_psi0) THEN
1030 DEALLOCATE (dcdr_env%apt_el_dcdr_per_center)
1031 DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset)
1032 END IF
1033
1034 ! Full matrices
1035 CALL cp_fm_release(dcdr_env%dCR)
1036 CALL cp_fm_release(dcdr_env%dCR_prime)
1037 CALL cp_fm_release(dcdr_env%mo_coeff)
1038 CALL cp_fm_release(dcdr_env%chc)
1039 CALL cp_fm_release(dcdr_env%op_dR)
1040
1041 ! DBCSR matrices
1042 CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
1043 CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)
1044 CALL dbcsr_deallocate_matrix_set(dcdr_env%moments)
1045 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip)
1046 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1)
1047 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp)
1048 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp2)
1049 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s)
1050 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t)
1051 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1)
1052 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1)
1053 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant)
1054 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR)
1055 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis)
1056 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc)
1057 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1)
1058
1059 END SUBROUTINE dcdr_env_cleanup
1060
1061END MODULE qs_dcdr_utils
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:89
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
cp_result_methods::get_results
Definition cp_result_methods.F:43
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:744
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:494
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:768
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:901
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:704
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::current_orb_center_wannier
integer, parameter, public current_orb_center_wannier
Definition input_constants.F:501
input_constants::use_mom_ref_user
integer, parameter, public use_mom_ref_user
Definition input_constants.F:496
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition moments_utils.F:61
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_dcdr_utils
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition qs_dcdr_utils.F:13
qs_dcdr_utils::dcdr_env_cleanup
subroutine, public dcdr_env_cleanup(qs_env, dcdr_env)
Deallocate the dcdr environment.
Definition qs_dcdr_utils.F:990
qs_dcdr_utils::multiply_localization
subroutine, public multiply_localization(ao_matrix, mo_coeff, work, nmo, icenter, res)
Multiply (ao_matrix @ mo_coeff) and store the column icenter in res.
Definition qs_dcdr_utils.F:113
qs_dcdr_utils::dcdr_print
subroutine, public dcdr_print(dcdr_env, qs_env)
Print the APT and sum rules.
Definition qs_dcdr_utils.F:397
qs_dcdr_utils::dcdr_read_restart
subroutine, public dcdr_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_read_restart.
Definition qs_dcdr_utils.F:148
qs_dcdr_utils::shift_wannier_into_cell
subroutine, public shift_wannier_into_cell(r, cell, r_shifted)
...
Definition qs_dcdr_utils.F:526
qs_dcdr_utils::get_loc_setting
subroutine, public get_loc_setting(dcdr_env, qs_env)
...
Definition qs_dcdr_utils.F:553
qs_dcdr_utils::dcdr_write_restart
subroutine, public dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
Copied from linres_write_restart.
Definition qs_dcdr_utils.F:298
qs_dcdr_utils::dcdr_env_init
subroutine, public dcdr_env_init(dcdr_env, qs_env)
Initialize the dcdr environment.
Definition qs_dcdr_utils.F:652
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_kinetic
Calculation of kinetic energy matrix and forces.
Definition qs_kinetic.F:15
qs_kinetic::build_kinetic_matrix
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition qs_kinetic.F:101
qs_ks_types
Definition qs_ks_types.F:15
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_loc_types
New version of the module for the localization of the molecular orbitals This should be able to use d...
Definition qs_loc_types.F:25
qs_loc_types::get_qs_loc_env
subroutine, public get_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
Definition qs_loc_types.F:317
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:558
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::xstring
elemental subroutine, public xstring(string, ia, ib)
...
Definition string_utilities.F:3381
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_linres_types::dcdr_env_type
Definition qs_linres_types.F:247
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_loc_types::localized_wfn_control_type
A type that holds controlling information for the calculation of the spread of wfn and the optimizati...
Definition qs_loc_types.F:121
qs_loc_types::qs_loc_env_type
contains all the info needed by quickstep to calculate the spread of a selected set of orbitals and i...
Definition qs_loc_types.F:80
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67