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