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