(git:ddb311d)
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 vecbuffer = 0.0_dp
232 !
233 ! read headers
234 IF (rst_unit > 0) READ (rst_unit, iostat=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
235 CALL para_env%bcast(iostat)
236 CALL para_env%bcast(beta_tmp)
237 CALL para_env%bcast(lambda_tmp)
238 CALL para_env%bcast(nspins_tmp)
239 CALL para_env%bcast(nao_tmp)
240
241 ! check that the number nao, nmo and nspins are
242 ! the same as in the current mos
243 IF (nspins_tmp /= nspins) THEN
244 cpabort("nspins not consistent")
245 END IF
246 IF (nao_tmp /= nao) cpabort("nao not consistent")
247 ! check that it's the right file
248 ! the same as in the current mos
249 IF (lambda_tmp /= lambda) cpabort("lambda not consistent")
250 IF (beta_tmp /= beta) cpabort("beta not consistent")
251 !
252 DO ispin = 1, nspins
253 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
254 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
255 !
256 IF (rst_unit > 0) READ (rst_unit) nmo_tmp
257 CALL para_env%bcast(nmo_tmp)
258 IF (nmo_tmp /= nmo) cpabort("nmo not consistent")
259 !
260 ! read the response
261 DO i = 1, nmo, max(max_block, 1)
262 i_block = min(max_block, nmo - i + 1)
263 DO j = 1, i_block
264 IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
265 END DO
266 CALL para_env%bcast(vecbuffer)
267 CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
268 END DO
269 END DO
270
271 IF (iostat /= 0) THEN
272 IF (iounit > 0) WRITE (iounit, "(T2,A)") &
273 "LINRES| Restart file <"//trim(adjustl(filename))//"> not found"
274 END IF
275
276 DEALLOCATE (vecbuffer)
277
278 END IF
279
280 IF (para_env%is_source()) THEN
281 IF (file_exists) CALL close_file(unit_number=rst_unit)
282 END IF
283
284 CALL timestop(handle)
285
286 END SUBROUTINE dcdr_read_restart
287
288! **************************************************************************************************
289!> \brief Copied from linres_write_restart
290!> \param qs_env ...
291!> \param linres_section ...
292!> \param vec ...
293!> \param lambda ...
294!> \param beta ...
295!> \param tag ...
296!> \note Adapted from linres_read_restart (ED)
297!> Would be nice not to crash but to start from zero if the present file doesn't match.
298! **************************************************************************************************
299 SUBROUTINE dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
300 TYPE(qs_environment_type), POINTER :: qs_env
301 TYPE(section_vals_type), POINTER :: linres_section
302 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: vec
303 INTEGER, INTENT(IN) :: lambda, beta
304 CHARACTER(LEN=*) :: tag
305
306 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_write_restart'
307
308 CHARACTER(LEN=default_path_length) :: filename
309 CHARACTER(LEN=default_string_length) :: my_middle, my_pos, my_status
310 INTEGER :: handle, i, i_block, ia, ie, iounit, &
311 ispin, j, max_block, nao, nmo, nspins, &
312 rst_unit
313 REAL(kind=dp), DIMENSION(:, :), POINTER :: vecbuffer
314 TYPE(cp_fm_type), POINTER :: mo_coeff
315 TYPE(cp_logger_type), POINTER :: logger
316 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
317 TYPE(mp_para_env_type), POINTER :: para_env
318 TYPE(section_vals_type), POINTER :: print_key
319
320 NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
321
322 CALL timeset(routinen, handle)
323
324 logger => cp_get_default_logger()
325
326 IF (btest(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
327 used_print_key=print_key), &
328 cp_p_file)) THEN
329
330 iounit = cp_print_key_unit_nr(logger, linres_section, &
331 "PRINT%PROGRAM_RUN_INFO", extension=".Log")
332
333 CALL get_qs_env(qs_env=qs_env, &
334 mos=mos, &
335 para_env=para_env)
336
337 nspins = SIZE(mos)
338
339 my_status = "REPLACE"
340 my_pos = "REWIND"
341 CALL xstring(tag, ia, ie)
342 my_middle = "RESTART-"//tag(ia:ie)//trim("-")//trim(adjustl(cp_to_string(beta))) &
343 //trim("-")//trim(adjustl(cp_to_string(lambda)))
344 rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
345 extension=".lr", middle_name=trim(my_middle), file_status=trim(my_status), &
346 file_position=trim(my_pos), file_action="WRITE", file_form="UNFORMATTED")
347
348 filename = cp_print_key_generate_filename(logger, print_key, &
349 extension=".lr", middle_name=trim(my_middle), my_local=.false.)
350
351 IF (iounit > 0) THEN
352 WRITE (unit=iounit, fmt="(T2,A)") &
353 "LINRES| Writing response functions to the restart file <"//trim(adjustl(filename))//">"
354 END IF
355
356 !
357 ! write data to file
358 ! use the scalapack block size as a default for buffering columns
359 CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
360 CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
361 ALLOCATE (vecbuffer(nao, max_block))
362
363 IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
364
365 DO ispin = 1, nspins
366 CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
367
368 IF (rst_unit > 0) WRITE (rst_unit) nmo
369
370 DO i = 1, nmo, max(max_block, 1)
371 i_block = min(max_block, nmo - i + 1)
372 CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
373 ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
374 ! to old ones, and in cases where max_block is different between runs, as might happen during
375 ! restarts with a different number of CPUs
376 DO j = 1, i_block
377 IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
378 END DO
379 END DO
380 END DO
381
382 DEALLOCATE (vecbuffer)
383
384 CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
385 "PRINT%RESTART")
386 END IF
387
388 CALL timestop(handle)
389
390 END SUBROUTINE dcdr_write_restart
391
392! **************************************************************************************************
393!> \brief Print the APT and sum rules
394!> \param dcdr_env ...
395!> \param qs_env ...
396!> \author Edward Ditler
397! **************************************************************************************************
398 SUBROUTINE dcdr_print(dcdr_env, qs_env)
399 TYPE(dcdr_env_type) :: dcdr_env
400 TYPE(qs_environment_type), POINTER :: qs_env
401
402 CHARACTER(LEN=default_string_length) :: description
403 INTEGER :: alpha, beta, delta, gamma, i, k, l, &
404 lambda, natom, nsubset, output_unit
405 REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
406 REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center_dcdr, apt_subset_dcdr
407 REAL(kind=dp), DIMENSION(3, 3) :: sum_rule_0, sum_rule_1, sum_rule_2
408 TYPE(cp_logger_type), POINTER :: logger
409 TYPE(cp_result_type), POINTER :: results
410 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
411 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
412 TYPE(section_vals_type), POINTER :: dcdr_section
413
414 NULLIFY (logger)
415
416 logger => cp_get_default_logger()
417 output_unit = cp_logger_get_default_io_unit(logger)
418
419 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
420
421 NULLIFY (particle_set)
422 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
423 natom = SIZE(particle_set)
424 nsubset = SIZE(molecule_set)
425
426 apt_el_dcdr => dcdr_env%apt_el_dcdr
427 apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr
428 apt_total_dcdr => dcdr_env%apt_total_dcdr
429 apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset
430 apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center
431
432 IF (dcdr_env%localized_psi0) THEN
433 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the final apt matrix per atom per subset'
434 DO k = 1, natom
435 DO l = 1, nsubset
436 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, *) 'APT | Subset', l
437 DO i = 1, 3
438 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,F15.6,F15.6,F15.6)") &
439 'APT | apt_subset ', i, apt_subset_dcdr(i, :, k, l)
440 END DO
441 END DO
442 END DO
443 END IF
444
445 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") &
446 'APT | Write the final apt matrix per atom (Position perturbation)'
447 DO l = 1, natom
448 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") &
449 'APT | Atom', l, ' - GAPT ', &
450 (apt_total_dcdr(1, 1, l) &
451 + apt_total_dcdr(2, 2, l) &
452 + apt_total_dcdr(3, 3, l))/3._dp
453 DO i = 1, 3
454 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
455 END DO
456 END DO
457
458 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix'
459 DO i = 1, 3
460 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
461 "(A,F15.6,F15.6,F15.6)") "APT | ", sum(apt_total_dcdr(i, :, :), dim=2)
462 END DO
463 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | End Write the final apt matrix'
464
465 ! Get the dipole
466 CALL get_qs_env(qs_env, results=results)
467 description = "[DIPOLE]"
468 CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3))
469
470 ! Sum rules [for all alpha, beta]
471 sum_rule_0 = 0._dp
472 sum_rule_1 = 0._dp
473 sum_rule_2 = 0._dp
474
475 DO alpha = 1, 3
476 DO beta = 1, 3
477 ! 0: sum_lambda apt(alpha, beta, lambda)
478 DO lambda = 1, natom
479 sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
480 + apt_total_dcdr(alpha, beta, lambda)
481 END DO
482
483 ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
484 DO gamma = 1, 3
485 sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
486 + levi_civita(alpha, beta, gamma)*dcdr_env%dipole_pos(gamma)
487 END DO
488
489 ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
490 DO lambda = 1, natom
491 DO gamma = 1, 3
492 DO delta = 1, 3
493 sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
494 + levi_civita(beta, gamma, delta) &
495 *particle_set(lambda)%r(gamma) &
496 *apt_total_dcdr(delta, alpha, lambda)
497 END DO
498 END DO
499 END DO
500
501 END DO ! beta
502 END DO ! alpha
503
504 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules"
505 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") &
506 "APT |", "Total APT", "Dipole", "R * APT"
507 DO alpha = 1, 3
508 DO beta = 1, 3
509 IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
510 "(A,I3,I3,F15.6,F15.6,F15.6)") &
511 "APT | ", &
512 alpha, beta, &
513 sum_rule_0(alpha, beta), &
514 sum_rule_1(alpha, beta), &
515 sum_rule_2(alpha, beta)
516 END DO
517 END DO
518
519 END SUBROUTINE dcdr_print
520
521! **************************************************************************************************
522!> \brief ...
523!> \param r ...
524!> \param cell ...
525!> \param r_shifted ...
526! **************************************************************************************************
527 SUBROUTINE shift_wannier_into_cell(r, cell, r_shifted)
528 REAL(dp), DIMENSION(3), INTENT(in) :: r
529 TYPE(cell_type), INTENT(in), POINTER :: cell
530 REAL(dp), DIMENSION(3), INTENT(out) :: r_shifted
531
532 INTEGER :: i
533 REAL(kind=dp), DIMENSION(3) :: abc
534
535 ! Only orthorombic cell for now
536 CALL get_cell(cell, abc=abc)
537
538 DO i = 1, 3
539 IF (r(i) < 0._dp) THEN
540 r_shifted(i) = r(i) + abc(i)
541 ELSE IF (r(i) > abc(i)) THEN
542 r_shifted(i) = r(i) - abc(i)
543 ELSE
544 r_shifted(i) = r(i)
545 END IF
546 END DO
547 END SUBROUTINE shift_wannier_into_cell
548
549! **************************************************************************************************
550!> \brief ...
551!> \param dcdr_env ...
552!> \param qs_env ...
553! **************************************************************************************************
554 SUBROUTINE get_loc_setting(dcdr_env, qs_env)
555 TYPE(dcdr_env_type) :: dcdr_env
556 TYPE(qs_environment_type), POINTER :: qs_env
557
558 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_loc_setting'
559
560 INTEGER :: handle, is, ispin, istate, max_states, &
561 nmo, nmoloc, nstate, nstate_list(2)
562 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
563 REAL(dp), DIMENSION(:, :), POINTER :: center_array
564 TYPE(linres_control_type), POINTER :: linres_control
565 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
566 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
567 TYPE(qs_loc_env_type), POINTER :: qs_loc_env
568 TYPE(section_vals_type), POINTER :: dcdr_section
569
570 CALL timeset(routinen, handle)
571
572 CALL get_qs_env(qs_env=qs_env, &
573 linres_control=linres_control, &
574 mos=mos)
575
576 ! Some checks
577 max_states = 0
578 CALL get_mo_set(mo_set=mos(1), nmo=nmo)
579 max_states = max(max_states, nmo)
580
581 ! check that the number of localized states is equal to the number of states
582 nmoloc = SIZE(linres_control%qs_loc_env%localized_wfn_control%centers_set(1)%array, 2)
583 IF (nmoloc /= nmo) THEN
584 cpabort("The number of localized functions is not equal to the number of states.")
585 END IF
586
587 ! which center for the orbitals shall we use
588 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
589 CALL section_vals_val_get(dcdr_section, "ORBITAL_CENTER", i_val=dcdr_env%orb_center)
590 SELECT CASE (dcdr_env%orb_center)
592 dcdr_env%orb_center_name = "WANNIER"
593 CASE DEFAULT
594 cpabort(" ")
595 END SELECT
596
597 qs_loc_env => linres_control%qs_loc_env
598 CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
599
600 ALLOCATE (dcdr_env%centers_set(dcdr_env%nspins))
601 ALLOCATE (dcdr_env%center_list(dcdr_env%nspins))
602 ALLOCATE (state_list(max_states, dcdr_env%nspins))
603 state_list(:, :) = huge(0)
604 nstate_list(:) = huge(0)
605
606 ! Build the state_list
607 DO ispin = 1, dcdr_env%nspins
608 center_array => localized_wfn_control%centers_set(ispin)%array
609 nstate = 0
610 DO istate = 1, SIZE(center_array, 2)
611 nstate = nstate + 1
612 state_list(nstate, ispin) = istate
613 END DO
614 nstate_list(ispin) = nstate
615
616 ! clustering the states
617 nstate = nstate_list(ispin)
618 dcdr_env%nstates(ispin) = nstate
619
620 ALLOCATE (dcdr_env%center_list(ispin)%array(2, nstate + 1))
621 ALLOCATE (dcdr_env%centers_set(ispin)%array(3, nstate))
622 dcdr_env%center_list(ispin)%array(:, :) = huge(0)
623 dcdr_env%centers_set(ispin)%array(:, :) = huge(0.0_dp)
624
625 center_array => localized_wfn_control%centers_set(ispin)%array
626
627 ! point to the psi0 centers
628 SELECT CASE (dcdr_env%orb_center)
630 ! use the wannier center as -center-
631 dcdr_env%nbr_center(ispin) = nstate
632 DO is = 1, nstate
633 istate = state_list(is, 1)
634 dcdr_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
635 dcdr_env%center_list(ispin)%array(1, is) = is
636 dcdr_env%center_list(ispin)%array(2, is) = istate
637 END DO
638 dcdr_env%center_list(ispin)%array(1, nstate + 1) = nstate + 1
639
640 CASE DEFAULT
641 cpabort("Unknown orbital center...")
642 END SELECT
643 END DO
644
645 CALL timestop(handle)
646 END SUBROUTINE get_loc_setting
647
648! **************************************************************************************************
649!> \brief Initialize the dcdr environment
650!> \param dcdr_env ...
651!> \param qs_env ...
652! **************************************************************************************************
653 SUBROUTINE dcdr_env_init(dcdr_env, qs_env)
654 TYPE(dcdr_env_type) :: dcdr_env
655 TYPE(qs_environment_type), POINTER :: qs_env
656
657 CHARACTER(LEN=*), PARAMETER :: routinen = 'dcdr_env_init'
658
659 INTEGER :: handle, homo, i, isize, ispin, j, jg, &
660 n_rep, nao, natom, nmo, nspins, &
661 nsubset, output_unit, reference, &
662 unit_number
663 INTEGER, DIMENSION(:), POINTER :: tmplist
664 LOGICAL :: explicit
665 REAL(kind=dp), DIMENSION(:), POINTER :: ref_point
666 TYPE(cp_fm_type) :: buf
667 TYPE(cp_fm_type), POINTER :: mo_coeff
668 TYPE(cp_logger_type), POINTER :: logger
669 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
670 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
671 TYPE(dft_control_type), POINTER :: dft_control
672 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
673 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
674 TYPE(mp_para_env_type), POINTER :: para_env
675 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
676 POINTER :: sab_all, sab_orb
677 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
678 TYPE(qs_ks_env_type), POINTER :: ks_env
679 TYPE(section_vals_type), POINTER :: dcdr_section, loc_section, lr_section
680
681 CALL timeset(routinen, handle)
682
683 ! Set up the logger
684 NULLIFY (logger, loc_section, dcdr_section, lr_section)
685 logger => cp_get_default_logger()
686 loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE")
687 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
688 dcdr_env%output_unit = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", &
689 extension=".data", middle_name="dcdr", log_filename=.false., &
690 file_position="REWIND", file_status="REPLACE")
691
692 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
693 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
694 extension=".linresLog")
695 unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
696
697 IF (output_unit > 0) THEN
698 WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***"
699 END IF
700
701 NULLIFY (ks_env, dft_control, sab_orb, sab_all, particle_set, molecule_set, matrix_s, matrix_ks, mos, para_env)
702 CALL get_qs_env(qs_env=qs_env, &
703 ks_env=ks_env, &
704 dft_control=dft_control, &
705 sab_orb=sab_orb, &
706 sab_all=sab_all, &
707 particle_set=particle_set, &
708 molecule_set=molecule_set, &
709 matrix_s=matrix_s, &
710 matrix_ks=matrix_ks, &
711 mos=mos, &
712 para_env=para_env)
713
714 natom = SIZE(particle_set)
715 nsubset = SIZE(molecule_set)
716 nspins = dft_control%nspins
717 dcdr_env%nspins = dft_control%nspins
718
719 NULLIFY (dcdr_env%matrix_s)
720 CALL build_overlap_matrix(ks_env, matrix_s=dcdr_env%matrix_s, &
721 matrix_name="OVERLAP MATRIX", &
722 nderivative=1, &
723 basis_type_a="ORB", &
724 basis_type_b="ORB", &
725 sab_nl=sab_orb)
726
727 NULLIFY (dcdr_env%matrix_t)
728 NULLIFY (matrix_p)
729 CALL kinetic_energy_matrix(qs_env, matrix_t=dcdr_env%matrix_t, matrix_p=matrix_p, &
730 matrix_name="KINETIC ENERGY MATRIX", &
731 calculate_forces=.false., &
732 basis_type="ORB", &
733 sab_orb=sab_orb, &
734 nderivative=1, &
735 eps_filter=dft_control%qs_control%eps_filter_matrix)
736
737! CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, &
738! matrix_name="KINETIC ENERGY MATRIX", &
739! basis_type="ORB", &
740! sab_nl=sab_orb, nderivative=1, &
741! eps_filter=dft_control%qs_control%eps_filter_matrix)
742
743 ! Get inputs
744 CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
745 CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
746 CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
747 CALL section_vals_val_get(dcdr_section, "Z_MATRIX_METHOD", l_val=dcdr_env%z_matrix_method)
748
749 dcdr_env%ref_point = 0._dp
750
751 ! List of atoms
752 NULLIFY (tmplist)
753 isize = 0
754 CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep)
755 IF (n_rep == 0) THEN
756 ALLOCATE (dcdr_env%list_of_atoms(natom))
757 DO jg = 1, natom
758 dcdr_env%list_of_atoms(jg) = jg
759 END DO
760 ELSE
761 DO jg = 1, n_rep
762 ALLOCATE (dcdr_env%list_of_atoms(isize))
763 CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist)
764 CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist))
765 dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
766 isize = SIZE(dcdr_env%list_of_atoms)
767 END DO
768 END IF
769
770 ! Reference point
771 IF (dcdr_env%localized_psi0) THEN
772 ! Get the Wannier localized wave functions and centers
773 CALL get_loc_setting(dcdr_env, qs_env)
774 ELSE
775 ! Get the reference point from the input
776 CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference)
777 CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit)
778 IF (explicit) THEN
779 CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point)
780 ELSE
781 IF (reference == use_mom_ref_user) &
782 cpabort("User-defined reference point should be given explicitly")
783 END IF
784
785 CALL get_reference_point(rpoint=dcdr_env%ref_point, qs_env=qs_env, &
786 reference=reference, &
787 ref_point=ref_point)
788 END IF
789
790 ! Helper matrix structs
791 NULLIFY (dcdr_env%aoao_fm_struct, &
792 dcdr_env%momo_fm_struct, &
793 dcdr_env%likemos_fm_struct, &
794 dcdr_env%homohomo_fm_struct)
795 CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, &
796 nao=nao, nmo=nmo, homo=homo)
797 CALL cp_fm_struct_create(dcdr_env%aoao_fm_struct, nrow_global=nao, &
798 ncol_global=nao, para_env=para_env, &
799 context=mo_coeff%matrix_struct%context)
800 CALL cp_fm_struct_create(dcdr_env%homohomo_fm_struct, nrow_global=homo, &
801 ncol_global=homo, para_env=para_env, &
802 context=mo_coeff%matrix_struct%context)
803 dcdr_env%nao = nao
804
805 ALLOCATE (dcdr_env%nmo(nspins))
806 ALLOCATE (dcdr_env%momo_fm_struct(nspins))
807 ALLOCATE (dcdr_env%likemos_fm_struct(nspins))
808 DO ispin = 1, nspins
809 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
810 nao=nao, nmo=nmo, homo=homo)
811 CALL cp_fm_struct_create(dcdr_env%momo_fm_struct(ispin)%struct, nrow_global=nmo, &
812 ncol_global=nmo, para_env=para_env, &
813 context=mo_coeff%matrix_struct%context)
814 CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct(ispin)%struct, &
815 template_fmstruct=mo_coeff%matrix_struct)
816 dcdr_env%nmo(ispin) = nmo
817 END DO
818
819 ! Fields of reals
820 ALLOCATE (dcdr_env%deltaR(3, natom))
821 ALLOCATE (dcdr_env%delta_basis_function(3, natom))
822 ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom))
823 ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom))
824 ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom))
825
826 dcdr_env%apt_el_dcdr = 0._dp
827 dcdr_env%apt_nuc_dcdr = 0._dp
828 dcdr_env%apt_total_dcdr = 0._dp
829
830 dcdr_env%deltaR = 0.0_dp
831 dcdr_env%delta_basis_function = 0._dp
832
833 ! Localization
834 IF (dcdr_env%localized_psi0) THEN
835 ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1)))
836 ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset))
837 ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset))
838 dcdr_env%apt_el_dcdr_per_center = 0._dp
839 dcdr_env%apt_el_dcdr_per_subset = 0._dp
840 dcdr_env%apt_subset = 0.0_dp
841 END IF
842
843 ! Full matrices
844 ALLOCATE (dcdr_env%mo_coeff(nspins))
845 ALLOCATE (dcdr_env%dCR(nspins))
846 ALLOCATE (dcdr_env%dCR_prime(nspins))
847 ALLOCATE (dcdr_env%chc(nspins))
848 ALLOCATE (dcdr_env%op_dR(nspins))
849
850 DO ispin = 1, nspins
851 CALL cp_fm_create(dcdr_env%dCR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.true.)
852 CALL cp_fm_create(dcdr_env%dCR_prime(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.true.)
853 CALL cp_fm_create(dcdr_env%mo_coeff(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.true.)
854 CALL cp_fm_create(dcdr_env%chc(ispin), dcdr_env%momo_fm_struct(ispin)%struct, set_zero=.true.)
855 CALL cp_fm_create(dcdr_env%op_dR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.true.)
856
857 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
858 CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
859 END DO
860
861 IF (dcdr_env%z_matrix_method) THEN
862 ALLOCATE (dcdr_env%matrix_m_alpha(3, nspins))
863 DO i = 1, 3
864 DO ispin = 1, nspins
865 CALL cp_fm_create(dcdr_env%matrix_m_alpha(i, ispin), dcdr_env%likemos_fm_struct(1)%struct)
866 CALL cp_fm_set_all(dcdr_env%matrix_m_alpha(i, ispin), 0.0_dp)
867 END DO
868 END DO
869 END IF
870
871 ! DBCSR matrices
872 NULLIFY (dcdr_env%hamiltonian1)
873 NULLIFY (dcdr_env%moments)
874 NULLIFY (dcdr_env%matrix_difdip)
875 NULLIFY (dcdr_env%matrix_core_charge_1)
876 NULLIFY (dcdr_env%matrix_s1)
877 NULLIFY (dcdr_env%matrix_t1)
878 NULLIFY (dcdr_env%matrix_apply_op_constant)
879 NULLIFY (dcdr_env%matrix_d_vhxc_dR)
880 NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis)
881 NULLIFY (dcdr_env%matrix_hc)
882 NULLIFY (dcdr_env%matrix_ppnl_1)
883 CALL dbcsr_allocate_matrix_set(dcdr_env%perturbed_dm_correction, nspins)
884 CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins)
885 CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3)
886 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3)
887 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3)
888 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3)
889 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp2, 3)
890 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4)
891 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4)
892 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, nspins)
893 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3, nspins)
894 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6)
895 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3)
896 CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3)
897
898 ! temporary no_symmetry matrix:
899 DO i = 1, 3
900 CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp(i)%matrix)
901 CALL dbcsr_create(dcdr_env%matrix_nosym_temp(i)%matrix, template=matrix_ks(1)%matrix, &
902 matrix_type=dbcsr_type_no_symmetry)
903 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all)
904 CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
905
906 CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp2(i)%matrix)
907 CALL dbcsr_create(dcdr_env%matrix_nosym_temp2(i)%matrix, template=matrix_ks(1)%matrix, &
908 matrix_type=dbcsr_type_no_symmetry)
909 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp2(i)%matrix, sab_all)
910 CALL dbcsr_set(dcdr_env%matrix_nosym_temp2(i)%matrix, 0._dp)
911 END DO
912
913 ! moments carry the result of build_local_moment_matrix
914 DO i = 1, 3
915 CALL dbcsr_init_p(dcdr_env%moments(i)%matrix)
916 CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments")
917 CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp)
918 END DO
919 CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp])
920
921 DO i = 1, 3
922 DO j = 1, 3
923 CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix)
924 CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
925 CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp)
926 END DO
927 END DO
928
929 DO ispin = 1, nspins
930 CALL dbcsr_init_p(dcdr_env%hamiltonian1(ispin)%matrix)
931 CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction(ispin)%matrix)
932 CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(ispin)%matrix)
933 CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(ispin)%matrix, matrix_ks(1)%matrix)
934 CALL dbcsr_copy(dcdr_env%perturbed_dm_correction(ispin)%matrix, matrix_ks(1)%matrix)
935 END DO
936
937 ! overlap/kinetic matrix: s(1) normal overlap matrix;
938 ! s(2:4) derivatives wrt. nuclear coordinates
939 CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix)
940 CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix)
941
942 CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix)
943 CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix)
944
945 DO i = 2, 4
946 CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix)
947 CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
948
949 CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix)
950 CALL dbcsr_copy(dcdr_env%matrix_t1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
951 END DO
952
953 ! j=1...3: derivative wrt nucleus A, 4...6: wrt nucleus B
954 DO ispin = 1, nspins
955 DO j = 1, 6
956 CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix)
957 CALL dbcsr_copy(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix, dcdr_env%matrix_s1(1)%matrix)
958 END DO
959 END DO
960
961 DO i = 1, 3
962 CALL dbcsr_init_p(dcdr_env%matrix_hc(i)%matrix)
963 CALL dbcsr_create(dcdr_env%matrix_hc(i)%matrix, template=matrix_ks(1)%matrix, &
964 matrix_type=dbcsr_type_symmetric)
965 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb)
966 CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
967 END DO
968
969 DO i = 1, 3
970 CALL dbcsr_init_p(dcdr_env%matrix_ppnl_1(i)%matrix)
971 CALL dbcsr_create(dcdr_env%matrix_ppnl_1(i)%matrix, template=matrix_ks(1)%matrix, &
972 matrix_type=dbcsr_type_symmetric)
973 CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb)
974 CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
975 END DO
976
977 DO i = 1, 3
978 DO ispin = 1, nspins
979 CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix)
980 CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, dcdr_env%matrix_s1(1)%matrix)
981 END DO
982
983 CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix)
984 CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix)
985 CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
986 END DO
987
988 ! CHC
989 DO ispin = 1, nspins
990 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
991 CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
992
993 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, nmo)
994 ! chc = mo * matrix_ks * mo
995 CALL parallel_gemm('T', 'N', nmo, nmo, nao, &
996 1.0_dp, mo_coeff, buf, &
997 0.0_dp, dcdr_env%chc(ispin))
998
999 CALL cp_fm_release(buf)
1000 END DO
1001
1002 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1003 "PRINT%PROGRAM_RUN_INFO")
1004
1005 CALL timestop(handle)
1006 END SUBROUTINE dcdr_env_init
1007
1008! **************************************************************************************************
1009!> \brief Deallocate the dcdr environment
1010!> \param qs_env ...
1011!> \param dcdr_env ...
1012! **************************************************************************************************
1013 SUBROUTINE dcdr_env_cleanup(qs_env, dcdr_env)
1014
1015 TYPE(qs_environment_type), POINTER :: qs_env
1016 TYPE(dcdr_env_type) :: dcdr_env
1017
1018 INTEGER :: ispin
1019 TYPE(cp_logger_type), POINTER :: logger
1020 TYPE(section_vals_type), POINTER :: dcdr_section
1021
1022 ! Destroy the logger
1023 logger => cp_get_default_logger()
1024 dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
1025 CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT")
1026
1027 DEALLOCATE (dcdr_env%list_of_atoms)
1028
1029 CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct)
1030 CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct)
1031 DO ispin = 1, dcdr_env%nspins
1032 CALL cp_fm_struct_release(dcdr_env%momo_fm_struct(ispin)%struct)
1033 CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct(ispin)%struct)
1034 END DO
1035 DEALLOCATE (dcdr_env%momo_fm_struct)
1036 DEALLOCATE (dcdr_env%likemos_fm_struct)
1037
1038 DEALLOCATE (dcdr_env%deltar)
1039 DEALLOCATE (dcdr_env%delta_basis_function)
1040
1041 IF (dcdr_env%localized_psi0) THEN
1042 ! DEALLOCATE (dcdr_env%psi0_order)
1043 DEALLOCATE (dcdr_env%centers_set(1)%array)
1044 DEALLOCATE (dcdr_env%center_list(1)%array)
1045 DEALLOCATE (dcdr_env%centers_set)
1046 DEALLOCATE (dcdr_env%center_list)
1047 DEALLOCATE (dcdr_env%apt_subset)
1048 END IF
1049
1050 DEALLOCATE (dcdr_env%apt_el_dcdr)
1051 DEALLOCATE (dcdr_env%apt_nuc_dcdr)
1052 DEALLOCATE (dcdr_env%apt_total_dcdr)
1053 IF (dcdr_env%localized_psi0) THEN
1054 DEALLOCATE (dcdr_env%apt_el_dcdr_per_center)
1055 DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset)
1056 END IF
1057
1058 ! Full matrices
1059 CALL cp_fm_release(dcdr_env%dCR)
1060 CALL cp_fm_release(dcdr_env%dCR_prime)
1061 CALL cp_fm_release(dcdr_env%mo_coeff)
1062 CALL cp_fm_release(dcdr_env%chc)
1063 CALL cp_fm_release(dcdr_env%op_dR)
1064
1065 IF (dcdr_env%z_matrix_method) THEN
1066 CALL cp_fm_release(dcdr_env%matrix_m_alpha)
1067 END IF
1068
1069 ! DBCSR matrices
1070 CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
1071 CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)
1072 CALL dbcsr_deallocate_matrix_set(dcdr_env%moments)
1073 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip)
1074 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1)
1075 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp)
1076 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp2)
1077 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s)
1078 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t)
1079 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1)
1080 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1)
1081 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant)
1082 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR)
1083 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis)
1084 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc)
1085 CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1)
1086
1087 END SUBROUTINE dcdr_env_cleanup
1088
1089END 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:200
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:311
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:122
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:504
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_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
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,...
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, sab_cneo, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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:60
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...