Line data Source code
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 :
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_core_matrices, ONLY: kinetic_energy_matrix
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type
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
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 864 : CALL timeset(routineN, handle)
124 :
125 : ! Multiply by the MO coefficients
126 864 : CALL cp_dbcsr_sm_fm_multiply(ao_matrix, mo_coeff, work, ncol=nmo)
127 :
128 : ! Only keep the icenter-th column
129 864 : CALL cp_fm_to_fm(work, res, 1, icenter, icenter)
130 :
131 : ! Reset the matrices
132 864 : CALL cp_fm_set_all(work, 0.0_dp)
133 :
134 864 : CALL timestop(handle)
135 864 : 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 18 : 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 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
163 : TYPE(cp_fm_type), POINTER :: mo_coeff
164 : TYPE(cp_logger_type), POINTER :: logger
165 18 : 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 18 : file_exists = .FALSE.
170 :
171 18 : CALL timeset(routineN, handle)
172 :
173 18 : NULLIFY (mos, para_env, logger, print_key, vecbuffer)
174 18 : logger => cp_get_default_logger()
175 :
176 : iounit = cp_print_key_unit_nr(logger, linres_section, &
177 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
178 :
179 : CALL get_qs_env(qs_env=qs_env, &
180 : para_env=para_env, &
181 18 : mos=mos)
182 :
183 18 : nspins = SIZE(mos)
184 :
185 18 : rst_unit = -1
186 18 : IF (para_env%is_source()) THEN
187 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
188 9 : n_rep_val=n_rep_val)
189 :
190 9 : CALL XSTRING(tag, ia, ie)
191 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
192 9 : //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
193 :
194 9 : IF (n_rep_val > 0) THEN
195 0 : CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
196 0 : CALL xstring(filename, ia, ie)
197 0 : 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 9 : print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
201 : filename = cp_print_key_generate_filename(logger, print_key, &
202 9 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
203 : END IF
204 9 : INQUIRE (FILE=filename, exist=file_exists)
205 : !
206 : ! open file
207 9 : 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 0 : unit_number=rst_unit)
214 :
215 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
216 0 : "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
217 : ELSE
218 9 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
219 9 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
220 : END IF
221 : END IF
222 :
223 18 : CALL para_env%bcast(file_exists)
224 :
225 18 : IF (file_exists) THEN
226 :
227 0 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
228 0 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
229 :
230 0 : ALLOCATE (vecbuffer(nao, max_block))
231 0 : vecbuffer = 0.0_dp
232 : !
233 : ! read headers
234 0 : IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
235 0 : CALL para_env%bcast(iostat)
236 0 : CALL para_env%bcast(beta_tmp)
237 0 : CALL para_env%bcast(lambda_tmp)
238 0 : CALL para_env%bcast(nspins_tmp)
239 0 : 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 0 : IF (nspins_tmp /= nspins) THEN
244 0 : CPABORT("nspins not consistent")
245 : END IF
246 0 : 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 0 : IF (lambda_tmp /= lambda) CPABORT("lambda not consistent")
250 0 : IF (beta_tmp /= beta) CPABORT("beta not consistent")
251 : !
252 0 : DO ispin = 1, nspins
253 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
254 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
255 : !
256 0 : IF (rst_unit > 0) READ (rst_unit) nmo_tmp
257 0 : CALL para_env%bcast(nmo_tmp)
258 0 : IF (nmo_tmp /= nmo) CPABORT("nmo not consistent")
259 : !
260 : ! read the response
261 0 : DO i = 1, nmo, MAX(max_block, 1)
262 0 : i_block = MIN(max_block, nmo - i + 1)
263 0 : DO j = 1, i_block
264 0 : IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
265 : END DO
266 0 : CALL para_env%bcast(vecbuffer)
267 0 : CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
268 : END DO
269 : END DO
270 :
271 0 : IF (iostat /= 0) THEN
272 0 : IF (iounit > 0) WRITE (iounit, "(T2,A)") &
273 0 : "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
274 : END IF
275 :
276 0 : DEALLOCATE (vecbuffer)
277 :
278 : END IF
279 :
280 18 : IF (para_env%is_source()) THEN
281 9 : IF (file_exists) CALL close_file(unit_number=rst_unit)
282 : END IF
283 :
284 18 : CALL timestop(handle)
285 :
286 18 : 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 18 : 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 18 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer
314 : TYPE(cp_fm_type), POINTER :: mo_coeff
315 : TYPE(cp_logger_type), POINTER :: logger
316 18 : 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 18 : NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
321 :
322 18 : CALL timeset(routineN, handle)
323 :
324 18 : logger => cp_get_default_logger()
325 :
326 18 : 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 18 : "PRINT%PROGRAM_RUN_INFO", extension=".Log")
332 :
333 : CALL get_qs_env(qs_env=qs_env, &
334 : mos=mos, &
335 18 : para_env=para_env)
336 :
337 18 : nspins = SIZE(mos)
338 :
339 18 : my_status = "REPLACE"
340 18 : my_pos = "REWIND"
341 18 : CALL XSTRING(tag, ia, ie)
342 : my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
343 18 : //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 18 : file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
347 :
348 : filename = cp_print_key_generate_filename(logger, print_key, &
349 18 : extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
350 :
351 18 : IF (iounit > 0) THEN
352 : WRITE (UNIT=iounit, FMT="(T2,A)") &
353 9 : "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 18 : CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
360 18 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
361 72 : ALLOCATE (vecbuffer(nao, max_block))
362 :
363 18 : IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
364 :
365 36 : DO ispin = 1, nspins
366 18 : CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
367 :
368 18 : IF (rst_unit > 0) WRITE (rst_unit) nmo
369 :
370 72 : DO i = 1, nmo, MAX(max_block, 1)
371 18 : i_block = MIN(max_block, nmo - i + 1)
372 18 : 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 108 : DO j = 1, i_block
377 90 : IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
378 : END DO
379 : END DO
380 : END DO
381 :
382 18 : DEALLOCATE (vecbuffer)
383 :
384 : CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
385 36 : "PRINT%RESTART")
386 : END IF
387 :
388 18 : CALL timestop(handle)
389 :
390 18 : 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 22 : 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 22 : REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
406 22 : 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 22 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
411 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
412 : TYPE(section_vals_type), POINTER :: dcdr_section
413 :
414 22 : NULLIFY (logger)
415 :
416 44 : logger => cp_get_default_logger()
417 22 : output_unit = cp_logger_get_default_io_unit(logger)
418 :
419 22 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
420 :
421 22 : NULLIFY (particle_set)
422 22 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
423 22 : natom = SIZE(particle_set)
424 22 : nsubset = SIZE(molecule_set)
425 :
426 22 : apt_el_dcdr => dcdr_env%apt_el_dcdr
427 22 : apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr
428 22 : apt_total_dcdr => dcdr_env%apt_total_dcdr
429 22 : apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset
430 22 : apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center
431 :
432 22 : IF (dcdr_env%localized_psi0) THEN
433 4 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the final apt matrix per atom per subset'
434 16 : DO k = 1, natom
435 52 : DO l = 1, nsubset
436 36 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, *) 'APT | Subset', l
437 156 : DO i = 1, 3
438 108 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,F15.6,F15.6,F15.6)") &
439 90 : 'APT | apt_subset ', i, apt_subset_dcdr(i, :, k, l)
440 : END DO
441 : END DO
442 : END DO
443 : END IF
444 :
445 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") &
446 11 : 'APT | Write the final apt matrix per atom (Position perturbation)'
447 88 : DO l = 1, natom
448 66 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") &
449 33 : 'APT | Atom', l, ' - GAPT ', &
450 : (apt_total_dcdr(1, 1, l) &
451 : + apt_total_dcdr(2, 2, l) &
452 66 : + apt_total_dcdr(3, 3, l))/3._dp
453 286 : DO i = 1, 3
454 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)
455 : END DO
456 : END DO
457 :
458 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix'
459 88 : DO i = 1, 3
460 66 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
461 451 : "(A,F15.6,F15.6,F15.6)") "APT | ", SUM(apt_total_dcdr(i, :, :), dim=2)
462 : END DO
463 22 : 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 22 : CALL get_qs_env(qs_env, results=results)
467 22 : description = "[DIPOLE]"
468 22 : CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3))
469 :
470 : ! Sum rules [for all alpha, beta]
471 22 : sum_rule_0 = 0._dp
472 22 : sum_rule_1 = 0._dp
473 22 : sum_rule_2 = 0._dp
474 :
475 88 : DO alpha = 1, 3
476 286 : DO beta = 1, 3
477 : ! 0: sum_lambda apt(alpha, beta, lambda)
478 792 : DO lambda = 1, natom
479 : sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
480 792 : + apt_total_dcdr(alpha, beta, lambda)
481 : END DO
482 :
483 : ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
484 792 : DO gamma = 1, 3
485 : sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
486 792 : + 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 858 : DO lambda = 1, natom
491 2574 : DO gamma = 1, 3
492 7722 : 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 7128 : *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 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules"
505 22 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") &
506 11 : "APT |", "Total APT", "Dipole", "R * APT"
507 88 : DO alpha = 1, 3
508 286 : DO beta = 1, 3
509 198 : IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
510 : "(A,I3,I3,F15.6,F15.6,F15.6)") &
511 99 : "APT | ", &
512 99 : alpha, beta, &
513 99 : sum_rule_0(alpha, beta), &
514 99 : sum_rule_1(alpha, beta), &
515 264 : sum_rule_2(alpha, beta)
516 : END DO
517 : END DO
518 :
519 22 : END SUBROUTINE dcdr_print
520 :
521 : ! **************************************************************************************************
522 : !> \brief ...
523 : !> \param r ...
524 : !> \param cell ...
525 : !> \param r_shifted ...
526 : ! **************************************************************************************************
527 144 : 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 144 : CALL get_cell(cell, abc=abc)
537 :
538 576 : DO i = 1, 3
539 576 : IF (r(i) < 0._dp) THEN
540 234 : r_shifted(i) = r(i) + abc(i)
541 198 : ELSE IF (r(i) > abc(i)) THEN
542 0 : r_shifted(i) = r(i) - abc(i)
543 : ELSE
544 198 : r_shifted(i) = r(i)
545 : END IF
546 : END DO
547 144 : END SUBROUTINE shift_wannier_into_cell
548 :
549 : ! **************************************************************************************************
550 : !> \brief ...
551 : !> \param dcdr_env ...
552 : !> \param qs_env ...
553 : ! **************************************************************************************************
554 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
563 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
571 :
572 : CALL get_qs_env(qs_env=qs_env, &
573 : linres_control=linres_control, &
574 4 : mos=mos)
575 :
576 : ! Some checks
577 4 : max_states = 0
578 4 : CALL get_mo_set(mo_set=mos(1), nmo=nmo)
579 4 : max_states = MAX(max_states, nmo)
580 :
581 : ! check that the number of localized states is equal to the number of states
582 4 : nmoloc = SIZE(linres_control%qs_loc_env%localized_wfn_control%centers_set(1)%array, 2)
583 4 : IF (nmoloc /= nmo) THEN
584 0 : 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 4 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
589 4 : CALL section_vals_val_get(dcdr_section, "ORBITAL_CENTER", i_val=dcdr_env%orb_center)
590 8 : SELECT CASE (dcdr_env%orb_center)
591 : CASE (current_orb_center_wannier)
592 4 : dcdr_env%orb_center_name = "WANNIER"
593 : CASE DEFAULT
594 4 : CPABORT(" ")
595 : END SELECT
596 :
597 4 : qs_loc_env => linres_control%qs_loc_env
598 4 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
599 :
600 16 : ALLOCATE (dcdr_env%centers_set(dcdr_env%nspins))
601 12 : ALLOCATE (dcdr_env%center_list(dcdr_env%nspins))
602 16 : ALLOCATE (state_list(max_states, dcdr_env%nspins))
603 24 : state_list(:, :) = HUGE(0)
604 12 : nstate_list(:) = HUGE(0)
605 :
606 : ! Build the state_list
607 8 : DO ispin = 1, dcdr_env%nspins
608 4 : center_array => localized_wfn_control%centers_set(ispin)%array
609 4 : nstate = 0
610 20 : DO istate = 1, SIZE(center_array, 2)
611 16 : nstate = nstate + 1
612 20 : state_list(nstate, ispin) = istate
613 : END DO
614 : nstate_list(ispin) = nstate
615 :
616 : ! clustering the states
617 4 : nstate = nstate_list(ispin)
618 4 : dcdr_env%nstates(ispin) = nstate
619 :
620 12 : ALLOCATE (dcdr_env%center_list(ispin)%array(2, nstate + 1))
621 12 : ALLOCATE (dcdr_env%centers_set(ispin)%array(3, nstate))
622 64 : dcdr_env%center_list(ispin)%array(:, :) = HUGE(0)
623 68 : dcdr_env%centers_set(ispin)%array(:, :) = HUGE(0.0_dp)
624 :
625 4 : center_array => localized_wfn_control%centers_set(ispin)%array
626 :
627 : ! point to the psi0 centers
628 4 : SELECT CASE (dcdr_env%orb_center)
629 : CASE (current_orb_center_wannier)
630 : ! use the wannier center as -center-
631 4 : dcdr_env%nbr_center(ispin) = nstate
632 20 : DO is = 1, nstate
633 16 : istate = state_list(is, 1)
634 64 : dcdr_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
635 16 : dcdr_env%center_list(ispin)%array(1, is) = is
636 20 : dcdr_env%center_list(ispin)%array(2, is) = istate
637 : END DO
638 4 : dcdr_env%center_list(ispin)%array(1, nstate + 1) = nstate + 1
639 :
640 : CASE DEFAULT
641 4 : CPABORT("Unknown orbital center...")
642 : END SELECT
643 : END DO
644 :
645 4 : CALL timestop(handle)
646 12 : END SUBROUTINE get_loc_setting
647 :
648 : ! **************************************************************************************************
649 : !> \brief Initialize the dcdr environment
650 : !> \param dcdr_env ...
651 : !> \param qs_env ...
652 : ! **************************************************************************************************
653 24 : 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 24 : INTEGER, DIMENSION(:), POINTER :: tmplist
664 : LOGICAL :: explicit
665 24 : 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 24 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
670 24 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
671 : TYPE(dft_control_type), POINTER :: dft_control
672 24 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
673 24 : 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 24 : POINTER :: sab_all, sab_orb
677 24 : 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 24 : CALL timeset(routineN, handle)
682 :
683 : ! Set up the logger
684 24 : NULLIFY (logger, loc_section, dcdr_section, lr_section)
685 24 : logger => cp_get_default_logger()
686 24 : loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE")
687 24 : 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 24 : file_position="REWIND", file_status="REPLACE")
691 :
692 24 : 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 24 : extension=".linresLog")
695 24 : unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
696 :
697 24 : IF (output_unit > 0) THEN
698 12 : WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***"
699 : END IF
700 :
701 24 : 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 24 : para_env=para_env)
713 :
714 24 : natom = SIZE(particle_set)
715 24 : nsubset = SIZE(molecule_set)
716 24 : nspins = dft_control%nspins
717 24 : dcdr_env%nspins = dft_control%nspins
718 :
719 24 : 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 24 : sab_nl=sab_orb)
726 :
727 24 : NULLIFY (dcdr_env%matrix_t)
728 24 : 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 24 : 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 24 : CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
745 24 : CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
746 24 : CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
747 24 : CALL section_vals_val_get(dcdr_section, "Z_MATRIX_METHOD", l_val=dcdr_env%z_matrix_method)
748 :
749 96 : dcdr_env%ref_point = 0._dp
750 :
751 : ! List of atoms
752 24 : NULLIFY (tmplist)
753 24 : isize = 0
754 24 : CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep)
755 24 : IF (n_rep == 0) THEN
756 72 : ALLOCATE (dcdr_env%list_of_atoms(natom))
757 96 : DO jg = 1, natom
758 96 : dcdr_env%list_of_atoms(jg) = jg
759 : END DO
760 : ELSE
761 0 : DO jg = 1, n_rep
762 0 : ALLOCATE (dcdr_env%list_of_atoms(isize))
763 0 : CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist)
764 0 : CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist))
765 0 : dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
766 0 : isize = SIZE(dcdr_env%list_of_atoms)
767 : END DO
768 : END IF
769 :
770 : ! Reference point
771 24 : IF (dcdr_env%localized_psi0) THEN
772 : ! Get the Wannier localized wave functions and centers
773 4 : CALL get_loc_setting(dcdr_env, qs_env)
774 : ELSE
775 : ! Get the reference point from the input
776 20 : CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference)
777 20 : CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit)
778 20 : IF (explicit) THEN
779 0 : CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point)
780 : ELSE
781 20 : IF (reference == use_mom_ref_user) &
782 0 : 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 20 : ref_point=ref_point)
788 : END IF
789 :
790 : ! Helper matrix structs
791 : NULLIFY (dcdr_env%aoao_fm_struct, &
792 24 : dcdr_env%momo_fm_struct, &
793 24 : dcdr_env%likemos_fm_struct, &
794 24 : dcdr_env%homohomo_fm_struct)
795 : CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, &
796 24 : 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 24 : 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 24 : context=mo_coeff%matrix_struct%context)
803 24 : dcdr_env%nao = nao
804 :
805 72 : ALLOCATE (dcdr_env%nmo(nspins))
806 100 : ALLOCATE (dcdr_env%momo_fm_struct(nspins))
807 76 : ALLOCATE (dcdr_env%likemos_fm_struct(nspins))
808 52 : DO ispin = 1, nspins
809 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
810 28 : 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 28 : context=mo_coeff%matrix_struct%context)
814 : CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct(ispin)%struct, &
815 28 : template_fmstruct=mo_coeff%matrix_struct)
816 80 : dcdr_env%nmo(ispin) = nmo
817 : END DO
818 :
819 : ! Fields of reals
820 72 : ALLOCATE (dcdr_env%deltaR(3, natom))
821 48 : ALLOCATE (dcdr_env%delta_basis_function(3, natom))
822 72 : ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom))
823 48 : ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom))
824 48 : ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom))
825 :
826 960 : dcdr_env%apt_el_dcdr = 0._dp
827 960 : dcdr_env%apt_nuc_dcdr = 0._dp
828 960 : dcdr_env%apt_total_dcdr = 0._dp
829 :
830 312 : dcdr_env%deltaR = 0.0_dp
831 312 : dcdr_env%delta_basis_function = 0._dp
832 :
833 : ! Localization
834 24 : IF (dcdr_env%localized_psi0) THEN
835 16 : ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1)))
836 16 : ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset))
837 12 : ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset))
838 644 : dcdr_env%apt_el_dcdr_per_center = 0._dp
839 484 : dcdr_env%apt_el_dcdr_per_subset = 0._dp
840 484 : dcdr_env%apt_subset = 0.0_dp
841 : END IF
842 :
843 : ! Full matrices
844 100 : ALLOCATE (dcdr_env%mo_coeff(nspins))
845 76 : ALLOCATE (dcdr_env%dCR(nspins))
846 76 : ALLOCATE (dcdr_env%dCR_prime(nspins))
847 76 : ALLOCATE (dcdr_env%chc(nspins))
848 76 : ALLOCATE (dcdr_env%op_dR(nspins))
849 :
850 52 : DO ispin = 1, nspins
851 28 : CALL cp_fm_create(dcdr_env%dCR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.TRUE.)
852 28 : CALL cp_fm_create(dcdr_env%dCR_prime(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.TRUE.)
853 28 : CALL cp_fm_create(dcdr_env%mo_coeff(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.TRUE.)
854 28 : CALL cp_fm_create(dcdr_env%chc(ispin), dcdr_env%momo_fm_struct(ispin)%struct, set_zero=.TRUE.)
855 28 : CALL cp_fm_create(dcdr_env%op_dR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct, set_zero=.TRUE.)
856 :
857 28 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
858 52 : CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
859 : END DO
860 :
861 24 : IF (dcdr_env%z_matrix_method) THEN
862 36 : ALLOCATE (dcdr_env%matrix_m_alpha(3, nspins))
863 16 : DO i = 1, 3
864 34 : DO ispin = 1, nspins
865 18 : CALL cp_fm_create(dcdr_env%matrix_m_alpha(i, ispin), dcdr_env%likemos_fm_struct(1)%struct)
866 30 : 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 24 : NULLIFY (dcdr_env%hamiltonian1)
873 24 : NULLIFY (dcdr_env%moments)
874 24 : NULLIFY (dcdr_env%matrix_difdip)
875 24 : NULLIFY (dcdr_env%matrix_core_charge_1)
876 24 : NULLIFY (dcdr_env%matrix_s1)
877 24 : NULLIFY (dcdr_env%matrix_t1)
878 24 : NULLIFY (dcdr_env%matrix_apply_op_constant)
879 24 : NULLIFY (dcdr_env%matrix_d_vhxc_dR)
880 24 : NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis)
881 24 : NULLIFY (dcdr_env%matrix_hc)
882 24 : NULLIFY (dcdr_env%matrix_ppnl_1)
883 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%perturbed_dm_correction, nspins)
884 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins)
885 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3)
886 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3)
887 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3)
888 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3)
889 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp2, 3)
890 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4)
891 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4)
892 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, nspins)
893 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3, nspins)
894 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6)
895 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3)
896 24 : CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3)
897 :
898 : ! temporary no_symmetry matrix:
899 96 : DO i = 1, 3
900 72 : 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 72 : matrix_type=dbcsr_type_no_symmetry)
903 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all)
904 72 : CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
905 :
906 72 : 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 72 : matrix_type=dbcsr_type_no_symmetry)
909 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp2(i)%matrix, sab_all)
910 96 : 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 96 : DO i = 1, 3
915 72 : CALL dbcsr_init_p(dcdr_env%moments(i)%matrix)
916 72 : CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments")
917 96 : CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp)
918 : END DO
919 24 : CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp])
920 :
921 96 : DO i = 1, 3
922 312 : DO j = 1, 3
923 216 : CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix)
924 216 : CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
925 288 : CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp)
926 : END DO
927 : END DO
928 :
929 52 : DO ispin = 1, nspins
930 28 : CALL dbcsr_init_p(dcdr_env%hamiltonian1(ispin)%matrix)
931 28 : CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction(ispin)%matrix)
932 28 : CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(ispin)%matrix)
933 28 : CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(ispin)%matrix, matrix_ks(1)%matrix)
934 52 : 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 24 : CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix)
940 24 : CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix)
941 :
942 24 : CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix)
943 24 : CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix)
944 :
945 96 : DO i = 2, 4
946 72 : CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix)
947 72 : CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
948 :
949 72 : CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix)
950 96 : 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 52 : DO ispin = 1, nspins
955 220 : DO j = 1, 6
956 168 : CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix)
957 196 : 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 96 : DO i = 1, 3
962 72 : 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 72 : matrix_type=dbcsr_type_symmetric)
965 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb)
966 96 : CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
967 : END DO
968 :
969 96 : DO i = 1, 3
970 72 : 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 72 : matrix_type=dbcsr_type_symmetric)
973 72 : CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb)
974 96 : CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
975 : END DO
976 :
977 96 : DO i = 1, 3
978 156 : DO ispin = 1, nspins
979 84 : CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix)
980 156 : CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, dcdr_env%matrix_s1(1)%matrix)
981 : END DO
982 :
983 72 : CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix)
984 72 : CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix)
985 96 : CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
986 : END DO
987 :
988 : ! CHC
989 52 : DO ispin = 1, nspins
990 28 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
991 28 : CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
992 :
993 28 : 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 28 : 0.0_dp, dcdr_env%chc(ispin))
998 :
999 80 : CALL cp_fm_release(buf)
1000 : END DO
1001 :
1002 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
1003 24 : "PRINT%PROGRAM_RUN_INFO")
1004 :
1005 24 : CALL timestop(handle)
1006 72 : END SUBROUTINE dcdr_env_init
1007 :
1008 : ! **************************************************************************************************
1009 : !> \brief Deallocate the dcdr environment
1010 : !> \param qs_env ...
1011 : !> \param dcdr_env ...
1012 : ! **************************************************************************************************
1013 24 : 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 24 : logger => cp_get_default_logger()
1024 24 : dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
1025 24 : CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT")
1026 :
1027 24 : DEALLOCATE (dcdr_env%list_of_atoms)
1028 :
1029 24 : CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct)
1030 24 : CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct)
1031 52 : DO ispin = 1, dcdr_env%nspins
1032 28 : CALL cp_fm_struct_release(dcdr_env%momo_fm_struct(ispin)%struct)
1033 52 : CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct(ispin)%struct)
1034 : END DO
1035 24 : DEALLOCATE (dcdr_env%momo_fm_struct)
1036 24 : DEALLOCATE (dcdr_env%likemos_fm_struct)
1037 :
1038 24 : DEALLOCATE (dcdr_env%deltar)
1039 24 : DEALLOCATE (dcdr_env%delta_basis_function)
1040 :
1041 24 : IF (dcdr_env%localized_psi0) THEN
1042 : ! DEALLOCATE (dcdr_env%psi0_order)
1043 4 : DEALLOCATE (dcdr_env%centers_set(1)%array)
1044 4 : DEALLOCATE (dcdr_env%center_list(1)%array)
1045 4 : DEALLOCATE (dcdr_env%centers_set)
1046 4 : DEALLOCATE (dcdr_env%center_list)
1047 4 : DEALLOCATE (dcdr_env%apt_subset)
1048 : END IF
1049 :
1050 24 : DEALLOCATE (dcdr_env%apt_el_dcdr)
1051 24 : DEALLOCATE (dcdr_env%apt_nuc_dcdr)
1052 24 : DEALLOCATE (dcdr_env%apt_total_dcdr)
1053 24 : IF (dcdr_env%localized_psi0) THEN
1054 4 : DEALLOCATE (dcdr_env%apt_el_dcdr_per_center)
1055 4 : DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset)
1056 : END IF
1057 :
1058 : ! Full matrices
1059 24 : CALL cp_fm_release(dcdr_env%dCR)
1060 24 : CALL cp_fm_release(dcdr_env%dCR_prime)
1061 24 : CALL cp_fm_release(dcdr_env%mo_coeff)
1062 24 : CALL cp_fm_release(dcdr_env%chc)
1063 24 : CALL cp_fm_release(dcdr_env%op_dR)
1064 :
1065 24 : IF (dcdr_env%z_matrix_method) THEN
1066 4 : CALL cp_fm_release(dcdr_env%matrix_m_alpha)
1067 : END IF
1068 :
1069 : ! DBCSR matrices
1070 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
1071 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)
1072 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%moments)
1073 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip)
1074 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1)
1075 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp)
1076 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp2)
1077 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s)
1078 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t)
1079 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1)
1080 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1)
1081 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant)
1082 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR)
1083 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis)
1084 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc)
1085 24 : CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1)
1086 :
1087 24 : END SUBROUTINE dcdr_env_cleanup
1088 :
1089 : END MODULE qs_dcdr_utils
|
