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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
19 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
20 : dbcsr_allocate_matrix_set,&
21 : dbcsr_deallocate_matrix_set
22 : USE cp_files, ONLY: close_file,&
23 : open_file
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_release
26 : USE cp_fm_types, ONLY: cp_fm_create,&
27 : cp_fm_get_info,&
28 : cp_fm_get_submatrix,&
29 : cp_fm_release,&
30 : cp_fm_set_all,&
31 : cp_fm_set_submatrix,&
32 : cp_fm_to_fm,&
33 : cp_fm_type
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_get_default_io_unit,&
36 : cp_logger_type,&
37 : cp_to_string
38 : USE cp_output_handling, ONLY: cp_p_file,&
39 : cp_print_key_finished_output,&
40 : cp_print_key_generate_filename,&
41 : cp_print_key_should_output,&
42 : cp_print_key_unit_nr
43 : USE cp_result_methods, ONLY: get_results
44 : USE cp_result_types, ONLY: cp_result_type
45 : USE dbcsr_api, ONLY: dbcsr_copy,&
46 : dbcsr_create,&
47 : dbcsr_init_p,&
48 : dbcsr_p_type,&
49 : dbcsr_set,&
50 : dbcsr_type,&
51 : dbcsr_type_no_symmetry,&
52 : dbcsr_type_symmetric
53 : USE input_constants, ONLY: current_orb_center_wannier,&
54 : use_mom_ref_user
55 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
56 : section_vals_type,&
57 : section_vals_val_get
58 : USE kinds, ONLY: default_path_length,&
59 : default_string_length,&
60 : dp
61 : USE memory_utilities, ONLY: reallocate
62 : USE message_passing, ONLY: mp_para_env_type
63 : USE molecule_types, ONLY: molecule_type
64 : USE moments_utils, ONLY: get_reference_point
65 : USE parallel_gemm_api, ONLY: parallel_gemm
66 : USE particle_types, ONLY: particle_type
67 : USE qs_environment_types, ONLY: get_qs_env,&
68 : qs_environment_type
69 : USE qs_kinetic, ONLY: build_kinetic_matrix
70 : USE qs_ks_types, ONLY: qs_ks_env_type
71 : USE qs_linres_types, ONLY: dcdr_env_type,&
72 : linres_control_type
73 : USE qs_loc_types, ONLY: get_qs_loc_env,&
74 : localized_wfn_control_type,&
75 : qs_loc_env_type
76 : USE qs_mo_types, ONLY: get_mo_set,&
77 : mo_set_type
78 : USE qs_moments, ONLY: build_local_moment_matrix
79 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
80 : USE qs_overlap, ONLY: build_overlap_matrix
81 : USE string_utilities, ONLY: xstring
82 : #include "./base/base_uses.f90"
83 :
84 : IMPLICIT NONE
85 :
86 : PRIVATE
87 : PUBLIC :: dcdr_env_cleanup, dcdr_env_init, dcdr_print, &
88 : get_loc_setting, shift_wannier_into_cell, &
89 : dcdr_write_restart, dcdr_read_restart, &
90 : multiply_localization
91 :
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_utils'
93 :
94 : REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = RESHAPE((/ &
95 : 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
96 : 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
97 : 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
98 : (/3, 3, 3/))
99 :
100 : CONTAINS
101 :
102 : ! **************************************************************************************************
103 : !> \brief Multiply (ao_matrix @ mo_coeff) and store the column icenter in res
104 : !> \param ao_matrix ...
105 : !> \param mo_coeff ...
106 : !> \param work Working space
107 : !> \param nmo ...
108 : !> \param icenter ...
109 : !> \param res ...
110 : !> \author Edward Ditler
111 : ! **************************************************************************************************
112 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 72 : DO i = 1, nmo, MAX(max_block, 1)
369 18 : i_block = MIN(max_block, nmo - i + 1)
370 18 : 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 108 : DO j = 1, i_block
375 90 : 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 10 : 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 10 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
404 10 : 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 10 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
409 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
410 : TYPE(section_vals_type), POINTER :: dcdr_section
411 :
412 10 : NULLIFY (logger)
413 :
414 20 : logger => cp_get_default_logger()
415 10 : output_unit = cp_logger_get_default_io_unit(logger)
416 :
417 10 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
418 :
419 10 : NULLIFY (particle_set)
420 10 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
421 10 : natom = SIZE(particle_set)
422 10 : nsubset = SIZE(molecule_set)
423 :
424 10 : apt_el_dcdr => dcdr_env%apt_el_dcdr
425 10 : apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr
426 10 : apt_total_dcdr => dcdr_env%apt_total_dcdr
427 10 : apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset
428 10 : apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center
429 :
430 10 : 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 10 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") &
444 5 : 'APT | Write the final apt matrix per atom (Position perturbation)'
445 40 : DO l = 1, natom
446 30 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") &
447 15 : 'APT | Atom', l, ' - GAPT ', &
448 : (apt_total_dcdr(1, 1, l) &
449 : + apt_total_dcdr(2, 2, l) &
450 30 : + apt_total_dcdr(3, 3, l))/3._dp
451 130 : DO i = 1, 3
452 120 : 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 10 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix'
457 40 : DO i = 1, 3
458 30 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
459 205 : "(A,F15.6,F15.6,F15.6)") "APT | ", SUM(apt_total_dcdr(i, :, :), dim=2)
460 : END DO
461 10 : 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 10 : CALL get_qs_env(qs_env, results=results)
465 10 : description = "[DIPOLE]"
466 10 : CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3))
467 :
468 : ! Sum rules [for all alpha, beta]
469 10 : sum_rule_0 = 0._dp
470 10 : sum_rule_1 = 0._dp
471 10 : sum_rule_2 = 0._dp
472 :
473 40 : DO alpha = 1, 3
474 130 : DO beta = 1, 3
475 : ! 0: sum_lambda apt(alpha, beta, lambda)
476 360 : DO lambda = 1, natom
477 : sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
478 360 : + apt_total_dcdr(alpha, beta, lambda)
479 : END DO
480 :
481 : ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
482 360 : DO gamma = 1, 3
483 : sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
484 360 : + 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 390 : DO lambda = 1, natom
489 1170 : DO gamma = 1, 3
490 3510 : 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 3240 : *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 10 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules"
503 10 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") &
504 5 : "APT |", "Total APT", "Dipole", "R * APT"
505 40 : DO alpha = 1, 3
506 130 : DO beta = 1, 3
507 90 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
508 : "(A,I3,I3,F15.6,F15.6,F15.6)") &
509 45 : "APT | ", &
510 45 : alpha, beta, &
511 45 : sum_rule_0(alpha, beta), &
512 45 : sum_rule_1(alpha, beta), &
513 120 : sum_rule_2(alpha, beta)
514 : END DO
515 : END DO
516 :
517 10 : 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 12 : 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 12 : INTEGER, DIMENSION(:), POINTER :: tmplist
662 : LOGICAL :: explicit
663 12 : 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 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
668 : TYPE(dft_control_type), POINTER :: dft_control
669 12 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
670 12 : 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 12 : POINTER :: sab_all, sab_orb
674 12 : 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 12 : CALL timeset(routineN, handle)
679 :
680 : ! Set up the logger
681 12 : NULLIFY (logger, loc_section, dcdr_section, lr_section)
682 12 : logger => cp_get_default_logger()
683 12 : loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE")
684 12 : 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 12 : file_position="REWIND", file_status="REPLACE")
688 :
689 12 : 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 12 : extension=".linresLog")
692 12 : unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
693 :
694 12 : IF (output_unit > 0) THEN
695 6 : WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***"
696 : END IF
697 :
698 12 : 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 12 : para_env=para_env)
710 :
711 12 : natom = SIZE(particle_set)
712 12 : nsubset = SIZE(molecule_set)
713 12 : nspins = dft_control%nspins
714 12 : dcdr_env%nspins = dft_control%nspins
715 :
716 12 : 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 12 : sab_nl=sab_orb)
723 :
724 12 : 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 12 : eps_filter=dft_control%qs_control%eps_filter_matrix)
730 :
731 : ! Get inputs
732 12 : CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
733 12 : CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
734 12 : CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
735 48 : dcdr_env%ref_point = 0._dp
736 :
737 : ! List of atoms
738 12 : NULLIFY (tmplist)
739 12 : isize = 0
740 12 : CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep)
741 12 : IF (n_rep == 0) THEN
742 36 : ALLOCATE (dcdr_env%list_of_atoms(natom))
743 48 : DO jg = 1, natom
744 48 : dcdr_env%list_of_atoms(jg) = jg
745 : END DO
746 : ELSE
747 0 : DO jg = 1, n_rep
748 0 : ALLOCATE (dcdr_env%list_of_atoms(isize))
749 0 : CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist)
750 0 : CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist))
751 0 : dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
752 0 : isize = SIZE(dcdr_env%list_of_atoms)
753 : END DO
754 : END IF
755 :
756 : ! Reference point
757 12 : IF (dcdr_env%localized_psi0) THEN
758 : ! Get the Wannier localized wave functions and centers
759 4 : CALL get_loc_setting(dcdr_env, qs_env)
760 : ELSE
761 : ! Get the reference point from the input
762 8 : CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference)
763 8 : CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit)
764 8 : IF (explicit) THEN
765 0 : CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point)
766 : ELSE
767 8 : IF (reference == use_mom_ref_user) &
768 0 : CPABORT("User-defined reference point should be given explicitly")
769 : END IF
770 :
771 : CALL get_reference_point(rpoint=dcdr_env%ref_point, qs_env=qs_env, &
772 : reference=reference, &
773 8 : ref_point=ref_point)
774 : END IF
775 :
776 : ! Helper matrix structs
777 : NULLIFY (dcdr_env%aoao_fm_struct, &
778 12 : dcdr_env%momo_fm_struct, &
779 12 : dcdr_env%likemos_fm_struct, &
780 12 : dcdr_env%homohomo_fm_struct)
781 : CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, &
782 12 : nao=nao, nmo=nmo, homo=homo)
783 : CALL cp_fm_struct_create(dcdr_env%aoao_fm_struct, nrow_global=nao, &
784 : ncol_global=nao, para_env=para_env, &
785 12 : context=mo_coeff%matrix_struct%context)
786 : CALL cp_fm_struct_create(dcdr_env%homohomo_fm_struct, nrow_global=homo, &
787 : ncol_global=homo, para_env=para_env, &
788 12 : context=mo_coeff%matrix_struct%context)
789 12 : dcdr_env%nao = nao
790 :
791 36 : ALLOCATE (dcdr_env%nmo(nspins))
792 50 : ALLOCATE (dcdr_env%momo_fm_struct(nspins))
793 38 : ALLOCATE (dcdr_env%likemos_fm_struct(nspins))
794 26 : DO ispin = 1, nspins
795 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
796 14 : nao=nao, nmo=nmo, homo=homo)
797 : CALL cp_fm_struct_create(dcdr_env%momo_fm_struct(ispin)%struct, nrow_global=nmo, &
798 : ncol_global=nmo, para_env=para_env, &
799 14 : context=mo_coeff%matrix_struct%context)
800 : CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct(ispin)%struct, &
801 14 : template_fmstruct=mo_coeff%matrix_struct)
802 40 : dcdr_env%nmo(ispin) = nmo
803 : END DO
804 :
805 : ! Fields of reals
806 36 : ALLOCATE (dcdr_env%deltaR(3, natom))
807 24 : ALLOCATE (dcdr_env%delta_basis_function(3, natom))
808 36 : ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom))
809 24 : ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom))
810 24 : ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom))
811 :
812 480 : dcdr_env%apt_el_dcdr = 0._dp
813 480 : dcdr_env%apt_nuc_dcdr = 0._dp
814 480 : dcdr_env%apt_total_dcdr = 0._dp
815 :
816 156 : dcdr_env%deltaR = 0.0_dp
817 156 : dcdr_env%delta_basis_function = 0._dp
818 :
819 : ! Localization
820 12 : IF (dcdr_env%localized_psi0) THEN
821 16 : ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1)))
822 16 : ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset))
823 12 : ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset))
824 644 : dcdr_env%apt_el_dcdr_per_center = 0._dp
825 484 : dcdr_env%apt_el_dcdr_per_subset = 0._dp
826 484 : dcdr_env%apt_subset = 0.0_dp
827 : END IF
828 :
829 : ! Full matrices
830 50 : ALLOCATE (dcdr_env%mo_coeff(nspins))
831 38 : ALLOCATE (dcdr_env%dCR(nspins))
832 38 : ALLOCATE (dcdr_env%dCR_prime(nspins))
833 38 : ALLOCATE (dcdr_env%chc(nspins))
834 38 : ALLOCATE (dcdr_env%op_dR(nspins))
835 :
836 26 : DO ispin = 1, nspins
837 14 : CALL cp_fm_create(dcdr_env%dCR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
838 14 : CALL cp_fm_create(dcdr_env%dCR_prime(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
839 14 : CALL cp_fm_create(dcdr_env%mo_coeff(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
840 14 : CALL cp_fm_create(dcdr_env%chc(ispin), dcdr_env%momo_fm_struct(ispin)%struct)
841 14 : CALL cp_fm_create(dcdr_env%op_dR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
842 :
843 14 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
844 26 : CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
845 : END DO
846 :
847 : ! DBCSR matrices
848 12 : NULLIFY (dcdr_env%hamiltonian1)
849 12 : NULLIFY (dcdr_env%moments)
850 12 : NULLIFY (dcdr_env%matrix_difdip)
851 12 : NULLIFY (dcdr_env%matrix_core_charge_1)
852 12 : NULLIFY (dcdr_env%matrix_s1)
853 12 : NULLIFY (dcdr_env%matrix_t1)
854 12 : NULLIFY (dcdr_env%matrix_apply_op_constant)
855 12 : NULLIFY (dcdr_env%matrix_d_vhxc_dR)
856 12 : NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis)
857 12 : NULLIFY (dcdr_env%matrix_hc)
858 12 : NULLIFY (dcdr_env%matrix_ppnl_1)
859 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%perturbed_dm_correction, nspins)
860 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins)
861 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3)
862 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3)
863 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3)
864 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3)
865 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp2, 3)
866 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4)
867 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4)
868 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, nspins)
869 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3, nspins)
870 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6)
871 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3)
872 12 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3)
873 :
874 : ! temporary no_symmetry matrix:
875 48 : DO i = 1, 3
876 36 : CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp(i)%matrix)
877 : CALL dbcsr_create(dcdr_env%matrix_nosym_temp(i)%matrix, template=matrix_ks(1)%matrix, &
878 36 : matrix_type=dbcsr_type_no_symmetry)
879 36 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all)
880 36 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
881 :
882 36 : CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp2(i)%matrix)
883 : CALL dbcsr_create(dcdr_env%matrix_nosym_temp2(i)%matrix, template=matrix_ks(1)%matrix, &
884 36 : matrix_type=dbcsr_type_no_symmetry)
885 36 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp2(i)%matrix, sab_all)
886 48 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp2(i)%matrix, 0._dp)
887 : END DO
888 :
889 : ! moments carry the result of build_local_moment_matrix
890 48 : DO i = 1, 3
891 36 : CALL dbcsr_init_p(dcdr_env%moments(i)%matrix)
892 36 : CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments")
893 48 : CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp)
894 : END DO
895 12 : CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp])
896 :
897 48 : DO i = 1, 3
898 156 : DO j = 1, 3
899 108 : CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix)
900 108 : CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
901 144 : CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp)
902 : END DO
903 : END DO
904 :
905 26 : DO ispin = 1, nspins
906 14 : CALL dbcsr_init_p(dcdr_env%hamiltonian1(ispin)%matrix)
907 14 : CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction(ispin)%matrix)
908 14 : CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(ispin)%matrix)
909 14 : CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(ispin)%matrix, matrix_ks(1)%matrix)
910 26 : CALL dbcsr_copy(dcdr_env%perturbed_dm_correction(ispin)%matrix, matrix_ks(1)%matrix)
911 : END DO
912 :
913 : ! overlap/kinetic matrix: s(1) normal overlap matrix;
914 : ! s(2:4) derivatives wrt. nuclear coordinates
915 12 : CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix)
916 12 : CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix)
917 :
918 12 : CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix)
919 12 : CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix)
920 :
921 48 : DO i = 2, 4
922 36 : CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix)
923 36 : CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
924 :
925 36 : CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix)
926 48 : CALL dbcsr_copy(dcdr_env%matrix_t1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
927 : END DO
928 :
929 : ! j=1...3: derivative wrt nucleus A, 4...6: wrt nucleus B
930 26 : DO ispin = 1, nspins
931 110 : DO j = 1, 6
932 84 : CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix)
933 98 : CALL dbcsr_copy(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix, dcdr_env%matrix_s1(1)%matrix)
934 : END DO
935 : END DO
936 :
937 48 : DO i = 1, 3
938 36 : CALL dbcsr_init_p(dcdr_env%matrix_hc(i)%matrix)
939 : CALL dbcsr_create(dcdr_env%matrix_hc(i)%matrix, template=matrix_ks(1)%matrix, &
940 36 : matrix_type=dbcsr_type_symmetric)
941 36 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb)
942 48 : CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
943 : END DO
944 :
945 48 : DO i = 1, 3
946 36 : CALL dbcsr_init_p(dcdr_env%matrix_ppnl_1(i)%matrix)
947 : CALL dbcsr_create(dcdr_env%matrix_ppnl_1(i)%matrix, template=matrix_ks(1)%matrix, &
948 36 : matrix_type=dbcsr_type_symmetric)
949 36 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb)
950 48 : CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
951 : END DO
952 :
953 48 : DO i = 1, 3
954 78 : DO ispin = 1, nspins
955 42 : CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix)
956 78 : CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, dcdr_env%matrix_s1(1)%matrix)
957 : END DO
958 :
959 36 : CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix)
960 36 : CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix)
961 48 : CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
962 : END DO
963 :
964 : ! CHC
965 26 : DO ispin = 1, nspins
966 14 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
967 14 : CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
968 :
969 14 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, nmo)
970 : ! chc = mo * matrix_ks * mo
971 : CALL parallel_gemm('T', 'N', nmo, nmo, nao, &
972 : 1.0_dp, mo_coeff, buf, &
973 14 : 0.0_dp, dcdr_env%chc(ispin))
974 :
975 40 : CALL cp_fm_release(buf)
976 : END DO
977 :
978 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
979 12 : "PRINT%PROGRAM_RUN_INFO")
980 :
981 12 : CALL timestop(handle)
982 36 : END SUBROUTINE dcdr_env_init
983 :
984 : ! **************************************************************************************************
985 : !> \brief Deallocate the dcdr environment
986 : !> \param qs_env ...
987 : !> \param dcdr_env ...
988 : ! **************************************************************************************************
989 12 : SUBROUTINE dcdr_env_cleanup(qs_env, dcdr_env)
990 :
991 : TYPE(qs_environment_type), POINTER :: qs_env
992 : TYPE(dcdr_env_type) :: dcdr_env
993 :
994 : INTEGER :: ispin
995 : TYPE(cp_logger_type), POINTER :: logger
996 : TYPE(section_vals_type), POINTER :: dcdr_section
997 :
998 : ! Destroy the logger
999 12 : logger => cp_get_default_logger()
1000 12 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
1001 12 : CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT")
1002 :
1003 12 : DEALLOCATE (dcdr_env%list_of_atoms)
1004 :
1005 12 : CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct)
1006 12 : CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct)
1007 26 : DO ispin = 1, dcdr_env%nspins
1008 14 : CALL cp_fm_struct_release(dcdr_env%momo_fm_struct(ispin)%struct)
1009 26 : CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct(ispin)%struct)
1010 : END DO
1011 12 : DEALLOCATE (dcdr_env%momo_fm_struct)
1012 12 : DEALLOCATE (dcdr_env%likemos_fm_struct)
1013 :
1014 12 : DEALLOCATE (dcdr_env%deltar)
1015 12 : DEALLOCATE (dcdr_env%delta_basis_function)
1016 :
1017 12 : IF (dcdr_env%localized_psi0) THEN
1018 : ! DEALLOCATE (dcdr_env%psi0_order)
1019 4 : DEALLOCATE (dcdr_env%centers_set(1)%array)
1020 4 : DEALLOCATE (dcdr_env%center_list(1)%array)
1021 4 : DEALLOCATE (dcdr_env%centers_set)
1022 4 : DEALLOCATE (dcdr_env%center_list)
1023 4 : DEALLOCATE (dcdr_env%apt_subset)
1024 : END IF
1025 :
1026 12 : DEALLOCATE (dcdr_env%apt_el_dcdr)
1027 12 : DEALLOCATE (dcdr_env%apt_nuc_dcdr)
1028 12 : DEALLOCATE (dcdr_env%apt_total_dcdr)
1029 12 : IF (dcdr_env%localized_psi0) THEN
1030 4 : DEALLOCATE (dcdr_env%apt_el_dcdr_per_center)
1031 4 : DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset)
1032 : END IF
1033 :
1034 : ! Full matrices
1035 12 : CALL cp_fm_release(dcdr_env%dCR)
1036 12 : CALL cp_fm_release(dcdr_env%dCR_prime)
1037 12 : CALL cp_fm_release(dcdr_env%mo_coeff)
1038 12 : CALL cp_fm_release(dcdr_env%chc)
1039 12 : CALL cp_fm_release(dcdr_env%op_dR)
1040 :
1041 : ! DBCSR matrices
1042 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
1043 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)
1044 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%moments)
1045 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip)
1046 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1)
1047 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp)
1048 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp2)
1049 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s)
1050 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t)
1051 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1)
1052 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1)
1053 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant)
1054 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR)
1055 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis)
1056 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc)
1057 12 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1)
1058 :
1059 12 : END SUBROUTINE dcdr_env_cleanup
1060 :
1061 : END MODULE qs_dcdr_utils
|
