(git:b279b6b)
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 
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_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 
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  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 
1061 END 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
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_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