Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Chemical shift calculation by dfpt
10 : !> Initialization of the nmr_env, creation of the special neighbor lists
11 : !> Perturbation Hamiltonians by application of the p and rxp oprtators to psi0
12 : !> Write output
13 : !> Deallocate everything
14 : !> \note
15 : !> The psi0 should be localized
16 : !> the Sebastiani method works within the assumption that the orbitals are
17 : !> completely contained in the simulation box
18 : !> \par History
19 : !> created 07-2005 [MI]
20 : !> \author MI
21 : ! **************************************************************************************************
22 : MODULE qs_linres_current_utils
23 : USE atomic_kind_types, ONLY: atomic_kind_type
24 : USE cell_types, ONLY: cell_type,&
25 : pbc
26 : USE cp_array_utils, ONLY: cp_2d_i_p_type,&
27 : cp_2d_r_p_type
28 : USE cp_control_types, ONLY: dft_control_type
29 : USE cp_dbcsr_api, ONLY: dbcsr_convert_offsets_to_sizes,&
30 : dbcsr_copy,&
31 : dbcsr_create,&
32 : dbcsr_distribution_type,&
33 : dbcsr_p_type,&
34 : dbcsr_set,&
35 : dbcsr_type_antisymmetric
36 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
37 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
38 : dbcsr_allocate_matrix_set,&
39 : dbcsr_deallocate_matrix_set
40 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
41 : cp_fm_scale_and_add
42 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43 : cp_fm_struct_release,&
44 : cp_fm_struct_type
45 : USE cp_fm_types, ONLY: cp_fm_create,&
46 : cp_fm_release,&
47 : cp_fm_set_all,&
48 : cp_fm_to_fm,&
49 : cp_fm_type
50 : USE cp_log_handling, ONLY: cp_get_default_logger,&
51 : cp_logger_type,&
52 : cp_to_string
53 : USE cp_output_handling, ONLY: cp_p_file,&
54 : cp_print_key_finished_output,&
55 : cp_print_key_should_output,&
56 : cp_print_key_unit_nr
57 : USE input_constants, ONLY: current_gauge_atom,&
58 : current_gauge_r,&
59 : current_gauge_r_and_step_func,&
60 : current_orb_center_atom,&
61 : current_orb_center_box,&
62 : current_orb_center_common,&
63 : current_orb_center_wannier,&
64 : ot_precond_full_all
65 : USE input_section_types, ONLY: section_get_lval,&
66 : section_vals_get_subs_vals,&
67 : section_vals_type,&
68 : section_vals_val_get
69 : USE kinds, ONLY: default_string_length,&
70 : dp
71 : USE memory_utilities, ONLY: reallocate
72 : USE message_passing, ONLY: mp_para_env_type
73 : USE particle_methods, ONLY: get_particle_set
74 : USE particle_types, ONLY: particle_type
75 : USE pw_env_types, ONLY: pw_env_get,&
76 : pw_env_type
77 : USE pw_methods, ONLY: pw_zero
78 : USE pw_pool_types, ONLY: pw_pool_type
79 : USE pw_types, ONLY: pw_c1d_gs_type,&
80 : pw_r3d_rs_type
81 : USE qs_environment_types, ONLY: get_qs_env,&
82 : qs_environment_type
83 : USE qs_kind_types, ONLY: qs_kind_type
84 : USE qs_linres_methods, ONLY: linres_read_restart,&
85 : linres_solver,&
86 : linres_write_restart
87 : USE qs_linres_op, ONLY: set_vecp
88 : USE qs_linres_types, ONLY: current_env_type,&
89 : deallocate_jrho_atom_set,&
90 : get_current_env,&
91 : init_jrho_atom_set,&
92 : jrho_atom_type,&
93 : linres_control_type,&
94 : set_current_env
95 : USE qs_loc_methods, ONLY: qs_print_cubes
96 : USE qs_loc_types, ONLY: get_qs_loc_env,&
97 : localized_wfn_control_type,&
98 : qs_loc_env_type
99 : USE qs_matrix_pools, ONLY: qs_matrix_pools_type
100 : USE qs_mo_types, ONLY: get_mo_set,&
101 : mo_set_type
102 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
103 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
104 : USE qs_p_env_types, ONLY: qs_p_env_type
105 : USE qs_rho_types, ONLY: qs_rho_clear,&
106 : qs_rho_create,&
107 : qs_rho_set
108 : USE realspace_grid_types, ONLY: rs_grid_release
109 : USE scf_control_types, ONLY: scf_control_type
110 : #include "./base/base_uses.f90"
111 :
112 : IMPLICIT NONE
113 :
114 : PRIVATE
115 : PUBLIC :: current_response, current_env_cleanup, current_env_init
116 :
117 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_current_utils'
118 :
119 : CONTAINS
120 :
121 : ! **************************************************************************************************
122 : !> \brief ...
123 : !> \param current_env ...
124 : !> \param p_env ...
125 : !> \param qs_env ...
126 : ! **************************************************************************************************
127 174 : SUBROUTINE current_response(current_env, p_env, qs_env)
128 : !
129 : TYPE(current_env_type) :: current_env
130 : TYPE(qs_p_env_type) :: p_env
131 : TYPE(qs_environment_type), POINTER :: qs_env
132 :
133 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_response'
134 :
135 : CHARACTER(LEN=default_string_length) :: my_pos
136 : INTEGER :: first_center, handle, i, icenter, idir, ii, iii, ispin, ist_true, istate, j, &
137 : jcenter, jstate, max_nbr_center, max_states, nao, natom, nbr_center(2), ncubes, nmo, &
138 : nspins, nstates(2), output_unit
139 174 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
140 174 : INTEGER, DIMENSION(:), POINTER :: list_cubes, row_blk_sizes
141 174 : INTEGER, DIMENSION(:, :, :), POINTER :: statetrueindex
142 : LOGICAL :: append_cube, should_stop
143 : REAL(dp) :: dk(3), dkl(3), dl(3)
144 174 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: dkl_vec_ii, dkl_vec_iii
145 174 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vecbuf_dklxp0
146 : TYPE(cell_type), POINTER :: cell
147 174 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
148 174 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
149 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
150 174 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_work_ii, fm_work_iii, h1_psi0, psi1
151 174 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: hpsi0_ptr, psi0_order
152 174 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: p_psi0, psi1_D, psi1_p, psi1_rxp, &
153 174 : rxp_psi0
154 : TYPE(cp_fm_type), POINTER :: mo_coeff
155 : TYPE(cp_logger_type), POINTER :: logger
156 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
157 174 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_p_ao
158 : TYPE(dft_control_type), POINTER :: dft_control
159 : TYPE(linres_control_type), POINTER :: linres_control
160 : TYPE(mp_para_env_type), POINTER :: para_env
161 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
162 174 : POINTER :: sab_orb
163 174 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
164 174 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
165 : TYPE(qs_matrix_pools_type), POINTER :: mpools
166 : TYPE(section_vals_type), POINTER :: current_section, lr_section, print_key
167 :
168 174 : CALL timeset(routineN, handle)
169 : !
170 174 : NULLIFY (cell, dft_control, linres_control, lr_section, current_section, &
171 174 : logger, mpools, mo_coeff, para_env, &
172 174 : tmp_fm_struct, &
173 174 : list_cubes, statetrueindex, centers_set, center_list, psi1_p, psi1_rxp, psi1_D, &
174 174 : p_psi0, rxp_psi0, psi0_order, op_p_ao, sab_orb, particle_set)
175 :
176 174 : logger => cp_get_default_logger()
177 174 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
178 : current_section => section_vals_get_subs_vals(qs_env%input, &
179 174 : "PROPERTIES%LINRES%CURRENT")
180 :
181 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
182 174 : extension=".linresLog")
183 174 : IF (output_unit > 0) THEN
184 : WRITE (UNIT=output_unit, FMT="(T10,A,/)") &
185 87 : "*** Self consistent optimization of the response wavefunctions ***"
186 : END IF
187 :
188 : CALL get_qs_env(qs_env=qs_env, &
189 : dft_control=dft_control, &
190 : mpools=mpools, cell=cell, &
191 : linres_control=linres_control, &
192 : sab_orb=sab_orb, &
193 : particle_set=particle_set, &
194 : qs_kind_set=qs_kind_set, &
195 : dbcsr_dist=dbcsr_dist, &
196 174 : para_env=para_env)
197 :
198 174 : nspins = dft_control%nspins
199 :
200 : CALL get_current_env(current_env=current_env, &
201 : nao=nao, &
202 : nstates=nstates, &
203 : centers_set=centers_set, &
204 : nbr_center=nbr_center, &
205 : center_list=center_list, &
206 : list_cubes=list_cubes, &
207 : statetrueindex=statetrueindex, &
208 : psi1_p=psi1_p, &
209 : psi1_rxp=psi1_rxp, &
210 : psi1_D=psi1_D, &
211 : p_psi0=p_psi0, &
212 : rxp_psi0=rxp_psi0, &
213 174 : psi0_order=psi0_order)
214 : !
215 : ! allocate the vectors
216 174 : IF (current_env%full) THEN
217 980 : ALLOCATE (psi1(nspins), h1_psi0(nspins))
218 348 : DO ispin = 1, nspins
219 206 : mo_coeff => psi0_order(ispin)
220 206 : NULLIFY (tmp_fm_struct)
221 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
222 : ncol_global=nstates(ispin), &
223 206 : context=mo_coeff%matrix_struct%context)
224 206 : CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
225 206 : CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
226 348 : CALL cp_fm_struct_release(tmp_fm_struct)
227 : END DO
228 : !
229 : ! prepare for allocation
230 142 : natom = SIZE(particle_set, 1)
231 426 : ALLOCATE (first_sgf(natom))
232 284 : ALLOCATE (last_sgf(natom))
233 : CALL get_particle_set(particle_set, qs_kind_set, &
234 : first_sgf=first_sgf, &
235 142 : last_sgf=last_sgf)
236 284 : ALLOCATE (row_blk_sizes(natom))
237 142 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
238 142 : DEALLOCATE (first_sgf, last_sgf)
239 : !
240 : ! rebuild the linear momentum matrices
241 142 : CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
242 142 : ALLOCATE (op_p_ao(1)%matrix)
243 : CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
244 : name="OP_P", &
245 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
246 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
247 142 : mutable_work=.TRUE.)
248 142 : CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
249 : !
250 : !
251 142 : DEALLOCATE (row_blk_sizes)
252 : !
253 : !
254 142 : CALL dbcsr_set(op_p_ao(1)%matrix, 0.0_dp)
255 426 : DO idir = 2, 3
256 284 : ALLOCATE (op_p_ao(idir)%matrix)
257 : CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
258 284 : "current_env%op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
259 426 : CALL dbcsr_set(op_p_ao(idir)%matrix, 0.0_dp)
260 : END DO
261 : !
262 142 : CALL build_lin_mom_matrix(qs_env, op_p_ao)
263 : !
264 : END IF
265 : !
266 : ! set response to zero
267 : !
268 696 : DO idir = 1, 3
269 1446 : DO ispin = 1, nspins
270 750 : CALL cp_fm_set_all(psi1_p(ispin, idir), 0.0_dp)
271 750 : CALL cp_fm_set_all(psi1_rxp(ispin, idir), 0.0_dp)
272 1272 : IF (current_env%full) CALL cp_fm_set_all(psi1_D(ispin, idir), 0.0_dp)
273 : END DO
274 : END DO
275 : !
276 : !
277 : !
278 : ! load restart file
279 : !
280 174 : first_center = 0
281 174 : IF (linres_control%linres_restart) THEN
282 24 : DO idir = 1, 3
283 : ! operator p
284 18 : CALL linres_read_restart(qs_env, lr_section, psi1_p(:, idir), idir, "nmr_p")
285 : ! operator rxp
286 24 : CALL linres_read_restart(qs_env, lr_section, psi1_rxp(:, idir), idir, "nmr_rxp")
287 : END DO
288 : END IF
289 : !
290 : !
291 : !
292 : !
293 174 : should_stop = .FALSE.
294 174 : current_env%all_pert_op_done = .FALSE.
295 : ! operator p
296 696 : DO idir = 1, 3
297 522 : IF (should_stop) EXIT
298 522 : IF (output_unit > 0) THEN
299 261 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator P_"//ACHAR(idir + 119)
300 : END IF
301 : !
302 : ! Initial guess for psi1
303 522 : hpsi0_ptr => p_psi0(:, idir)
304 : !
305 : !
306 522 : linres_control%converged = .FALSE.
307 522 : CALL linres_solver(p_env, qs_env, psi1_p(:, idir), hpsi0_ptr, psi0_order, output_unit, should_stop)
308 : !
309 : !
310 : ! print response functions
311 522 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
312 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
313 0 : ncubes = SIZE(list_cubes, 1)
314 0 : print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
315 0 : append_cube = section_get_lval(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%APPEND")
316 0 : my_pos = "REWIND"
317 0 : IF (append_cube) THEN
318 0 : my_pos = "APPEND"
319 : END IF
320 : !
321 :
322 0 : DO ispin = 1, nspins
323 : CALL qs_print_cubes(qs_env, psi1_p(ispin, idir), ncubes, list_cubes, &
324 : centers_set(ispin)%array, print_key, 'psi1_p', &
325 0 : idir=idir, ispin=ispin, file_position=my_pos)
326 : END DO ! ispin
327 : END IF ! print response functions
328 : !
329 : ! write restart file
330 696 : CALL linres_write_restart(qs_env, lr_section, psi1_p(:, idir), idir, "nmr_p")
331 : END DO ! idir
332 : !
333 : ! operator rxp
334 696 : DO idir = 1, 3
335 522 : IF (should_stop) EXIT
336 522 : IF (output_unit > 0) THEN
337 261 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator L_"//ACHAR(idir + 119)
338 : END IF
339 : !
340 : ! Initial guess for psi1
341 522 : hpsi0_ptr => rxp_psi0(:, idir)
342 : !
343 : !
344 522 : linres_control%converged = .FALSE.
345 522 : CALL linres_solver(p_env, qs_env, psi1_rxp(:, idir), hpsi0_ptr, psi0_order, output_unit, should_stop)
346 : !
347 : ! print response functions
348 522 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
349 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
350 0 : ncubes = SIZE(list_cubes, 1)
351 0 : print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
352 :
353 0 : DO ispin = 1, nspins
354 : CALL qs_print_cubes(qs_env, psi1_rxp(ispin, idir), ncubes, list_cubes, &
355 : centers_set(ispin)%array, print_key, 'psi1_rxp', &
356 0 : idir=idir, ispin=ispin, file_position=my_pos)
357 : END DO ! ispin
358 : END IF ! print response functions
359 : !
360 : ! write restart file
361 696 : CALL linres_write_restart(qs_env, lr_section, psi1_rxp(:, idir), idir, "nmr_rxp")
362 : END DO ! idir
363 174 : IF (.NOT. should_stop) current_env%all_pert_op_done = .TRUE.
364 : !
365 : ! operator D
366 174 : IF (current_env%full) THEN
367 142 : current_env%all_pert_op_done = .FALSE.
368 : !
369 : !
370 348 : DO ispin = 1, nspins
371 348 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
372 : END DO
373 : !
374 : ! The correction is state depedent a loop over the states is necessary
375 348 : max_nbr_center = MAXVAL(nbr_center(1:nspins))
376 348 : max_states = MAXVAL(nstates(1:nspins))
377 : !
378 0 : ALLOCATE (vecbuf_dklxp0(1, nao), fm_work_ii(nspins), fm_work_iii(nspins), &
379 1690 : dkl_vec_ii(max_states), dkl_vec_iii(max_states))
380 142 : vecbuf_dklxp0(1, nao) = 0.0_dp
381 : !
382 348 : DO ispin = 1, nspins
383 206 : nmo = nstates(ispin)
384 206 : mo_coeff => psi0_order(ispin)
385 206 : NULLIFY (tmp_fm_struct)
386 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
387 : ncol_global=nmo, para_env=para_env, &
388 206 : context=mo_coeff%matrix_struct%context)
389 :
390 206 : CALL cp_fm_create(fm_work_ii(ispin), tmp_fm_struct)
391 206 : CALL cp_fm_set_all(fm_work_ii(ispin), 0.0_dp)
392 206 : CALL cp_fm_create(fm_work_iii(ispin), tmp_fm_struct)
393 206 : CALL cp_fm_set_all(fm_work_iii(ispin), 0.0_dp)
394 348 : CALL cp_fm_struct_release(tmp_fm_struct)
395 : END DO ! ispin
396 : !
397 568 : DO idir = 1, 3
398 426 : IF (should_stop) EXIT
399 1044 : DO ispin = 1, nspins
400 1044 : CALL cp_fm_set_all(psi1_D(ispin, idir), 0.0_dp)
401 : END DO
402 426 : first_center = 0
403 426 : IF (linres_control%linres_restart) THEN
404 18 : CALL linres_read_restart(qs_env, lr_section, psi1_D(:, idir), idir, "nmr_dxp-", ind=first_center)
405 : END IF
406 426 : IF (first_center > 0) THEN
407 6 : IF (output_unit > 0) THEN
408 : WRITE (output_unit, "(T10,A,I4,A)")&
409 3 : & "Response to the perturbation operators (dk-dl)xp up to state ", &
410 6 : first_center, " have been read from restart"
411 : END IF
412 : END IF
413 : !
414 : ! here we run over the max number of states, then we need
415 : ! to be careful with overflow while doing uks calculations.
416 2886 : DO icenter = 1, max_nbr_center
417 : !
418 2460 : IF (should_stop) EXIT
419 2460 : IF (icenter > first_center) THEN
420 2436 : IF (output_unit > 0) THEN
421 : WRITE (output_unit, "(T10,A,I4,A)")&
422 1218 : & "Response to the perturbation operator (dk-dl)xp for -state- ", &
423 2436 : icenter, " in dir. "//ACHAR(idir + 119)
424 : END IF
425 : !
426 6306 : DO ispin = 1, nspins
427 3870 : nmo = nstates(ispin)
428 3870 : mo_coeff => psi0_order(ispin)
429 : !
430 : ! take care that no overflow can occur for uks
431 3870 : IF (icenter > nbr_center(ispin)) THEN
432 : !
433 : ! set h1_psi0 and psi1 to zero to avoid problems in linres_scf
434 276 : CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
435 276 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
436 276 : CYCLE
437 : END IF
438 : !
439 3594 : dkl_vec_ii(:) = 0.0_dp
440 3594 : dkl_vec_iii(:) = 0.0_dp
441 : !
442 3594 : ist_true = statetrueindex(idir, icenter, ispin)
443 : !
444 : ! the initial guess is the previous set of psi1, just optimized
445 3594 : CALL set_vecp(idir, ii, iii)
446 14376 : dk(1:3) = centers_set(ispin)%array(1:3, ist_true)
447 : !
448 27612 : DO jcenter = 1, nbr_center(ispin)
449 96072 : dl(1:3) = centers_set(ispin)%array(1:3, jcenter)
450 24018 : dkl = pbc(dl, dk, cell)
451 52962 : DO j = center_list(ispin)%array(1, jcenter), center_list(ispin)%array(1, jcenter + 1) - 1
452 25350 : jstate = center_list(ispin)%array(2, j)
453 25350 : dkl_vec_ii(jstate) = dkl(ii)
454 49368 : dkl_vec_iii(jstate) = dkl(iii)
455 : END DO
456 : END DO
457 : !
458 : ! First term
459 : ! Rescale the ground state orbitals by (dk-dl)_ii
460 3594 : CALL cp_fm_to_fm(mo_coeff, fm_work_ii(ispin))
461 3594 : CALL cp_fm_column_scale(fm_work_ii(ispin), dkl_vec_ii(1:nmo))
462 : !
463 : ! Apply the p_iii operator
464 : ! fm_work_iii = -p_iii * (dk-dl)_ii * C0
465 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(iii)%matrix, fm_work_ii(ispin), &
466 3594 : fm_work_iii(ispin), ncol=nmo, alpha=-1.0_dp)
467 : !
468 : ! Copy in h1_psi0
469 : ! h1_psi0_i = fm_work_iii
470 3594 : CALL cp_fm_to_fm(fm_work_iii(ispin), h1_psi0(ispin))
471 : !
472 : ! Second term
473 : ! Rescale the ground state orbitals by (dk-dl)_iii
474 3594 : CALL cp_fm_to_fm(mo_coeff, fm_work_iii(ispin))
475 3594 : CALL cp_fm_column_scale(fm_work_iii(ispin), dkl_vec_iii(1:nmo))
476 : !
477 : ! Apply the p_ii operator
478 : ! fm_work_ii = -p_ii * (dk-dl)_iii * C0
479 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(ii)%matrix, fm_work_iii(ispin), &
480 3594 : fm_work_ii(ispin), ncol=nmo, alpha=-1.0_dp)
481 : !
482 : ! Copy in h1_psi0
483 : ! h1_psi0_i = fm_work_iii - fm_work_ii
484 : CALL cp_fm_scale_and_add(1.0_dp, h1_psi0(ispin),&
485 9624 : & -1.0_dp, fm_work_ii(ispin))
486 : END DO
487 :
488 : !
489 : ! Optimize the response wavefunctions
490 2436 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
491 : !
492 2436 : IF (output_unit > 0) THEN
493 : WRITE (output_unit, "(T10,A,/)")&
494 1218 : & "Store the psi1 vector for the calculation of the response current density "
495 : END IF
496 : !
497 6306 : DO ispin = 1, nspins
498 : !
499 : ! take care that no overflow can occur for uks
500 3870 : IF (icenter > nbr_center(ispin)) CYCLE
501 : !
502 : ! need to reset those guys
503 3594 : ist_true = statetrueindex(idir, icenter, ispin)
504 7566 : DO i = center_list(ispin)%array(1, ist_true), center_list(ispin)%array(1, ist_true + 1) - 1
505 3972 : istate = center_list(ispin)%array(2, i)
506 : !
507 : ! the optimized wfns are copied in the fm
508 7566 : CALL cp_fm_to_fm(psi1(ispin), psi1_D(ispin, idir), 1, istate, istate)
509 : END DO
510 6306 : current_env%full_done(idir*ispin, icenter) = .TRUE.
511 : END DO ! ispin
512 : !
513 : ELSE
514 48 : DO ispin = 1, nspins
515 48 : current_env%full_done(idir*ispin, icenter) = .TRUE.
516 : END DO ! ispin
517 : END IF
518 2886 : CALL linres_write_restart(qs_env, lr_section, psi1_D(:, idir), idir, "nmr_dxp-", ind=icenter)
519 :
520 : END DO ! center
521 : !
522 : ! print response functions
523 426 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
524 142 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
525 0 : ncubes = SIZE(list_cubes, 1)
526 0 : print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
527 0 : DO ispin = 1, nspins
528 : CALL qs_print_cubes(qs_env, psi1_D(ispin, idir), &
529 : ncubes, list_cubes, centers_set(ispin)%array, print_key, 'psi1_D', &
530 0 : idir=idir, ispin=ispin, file_position=my_pos)
531 : END DO
532 : END IF ! print response functions
533 : !
534 : END DO ! idir
535 142 : IF (.NOT. should_stop) current_env%all_pert_op_done = .TRUE.
536 : !
537 : ! clean up
538 142 : CALL cp_fm_release(fm_work_ii)
539 142 : CALL cp_fm_release(fm_work_iii)
540 142 : DEALLOCATE (dkl_vec_ii, dkl_vec_iii, vecbuf_dklxp0)
541 :
542 : END IF
543 : !
544 : ! clean up
545 174 : IF (current_env%full) THEN
546 142 : CALL dbcsr_deallocate_matrix_set(op_p_ao)
547 142 : CALL cp_fm_release(psi1)
548 142 : CALL cp_fm_release(h1_psi0)
549 : END IF
550 : !
551 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
552 174 : & "PRINT%PROGRAM_RUN_INFO")
553 : !
554 174 : CALL timestop(handle)
555 : !
556 348 : END SUBROUTINE current_response
557 :
558 : ! **************************************************************************************************
559 :
560 : ! **************************************************************************************************
561 : !> \brief ...
562 : !> \param current_env ...
563 : !> \param qs_env ...
564 : ! **************************************************************************************************
565 174 : SUBROUTINE current_env_init(current_env, qs_env)
566 : !
567 : TYPE(current_env_type) :: current_env
568 : TYPE(qs_environment_type), POINTER :: qs_env
569 :
570 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_env_init'
571 :
572 : INTEGER :: handle, homo, i, iao, iatom, ibox, icenter, icount, idir, ii, ini, ir, is, ispin, &
573 : istate, istate2, istate_next, ix, iy, iz, j, jstate, k, max_nbr_center, max_states, n, &
574 : n_rep, nao, natom, nbr_box, ncubes, nmo, nspins, nstate, nstate_list(2), output_unit
575 174 : INTEGER, ALLOCATABLE, DIMENSION(:) :: buff, first_sgf, last_sgf
576 174 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
577 174 : INTEGER, DIMENSION(:), POINTER :: bounds, list, nbox, &
578 174 : selected_states_on_atom_list
579 : LOGICAL :: force_no_full, gapw, is0, &
580 : uniform_occupation
581 174 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: state_done
582 : REAL(dp) :: center(3), center2(3), dist, mdist, &
583 : r(3), rab(3)
584 174 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rbuff
585 174 : REAL(dp), DIMENSION(:), POINTER :: common_center
586 174 : REAL(dp), DIMENSION(:, :), POINTER :: center_array
587 174 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
588 : TYPE(cell_type), POINTER :: cell
589 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
590 : TYPE(cp_fm_type), POINTER :: mo_coeff
591 : TYPE(cp_logger_type), POINTER :: logger
592 : TYPE(dft_control_type), POINTER :: dft_control
593 174 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
594 : TYPE(linres_control_type), POINTER :: linres_control
595 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
596 174 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
597 : TYPE(mp_para_env_type), POINTER :: para_env
598 174 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
599 174 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
600 : TYPE(pw_env_type), POINTER :: pw_env
601 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
602 174 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
603 174 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
604 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
605 : TYPE(qs_matrix_pools_type), POINTER :: mpools
606 : TYPE(scf_control_type), POINTER :: scf_control
607 : TYPE(section_vals_type), POINTER :: current_section, lr_section
608 :
609 174 : CALL timeset(routineN, handle)
610 :
611 174 : NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, linres_control, scf_control, &
612 174 : logger, mos, mpools, current_section, particle_set, mo_coeff, &
613 174 : auxbas_pw_pool, pw_env, jrho1_atom_set, common_center, tmp_fm_struct, &
614 174 : para_env, qs_loc_env, localized_wfn_control, rho_g, rho_r)
615 : !
616 : !
617 : CALL get_qs_env(qs_env=qs_env, &
618 : atomic_kind_set=atomic_kind_set, &
619 : qs_kind_set=qs_kind_set, &
620 : cell=cell, &
621 : dft_control=dft_control, &
622 : linres_control=linres_control, &
623 : mos=mos, &
624 : mpools=mpools, &
625 : particle_set=particle_set, &
626 : pw_env=pw_env, &
627 : scf_control=scf_control, &
628 174 : para_env=para_env)
629 : !
630 174 : gapw = dft_control%qs_control%gapw
631 174 : nspins = dft_control%nspins
632 174 : natom = SIZE(particle_set, 1)
633 174 : CALL get_mo_set(mo_set=mos(1), nao=nao)
634 : !
635 174 : max_states = 0
636 424 : DO ispin = 1, nspins
637 250 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
638 424 : max_states = MAX(max_states, nmo)
639 : END DO
640 : !
641 : !
642 : !
643 : ! some checks
644 424 : DO ispin = 1, nspins
645 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, homo=homo, &
646 250 : uniform_occupation=uniform_occupation)
647 : !
648 : ! check that homo=nmo
649 250 : IF (nmo /= homo) CPABORT("nmo != homo")
650 : !
651 : ! check that the nbr of localized states is equal to the nbr of states
652 : ! nmoloc = SIZE(linres_control%localized_wfn_control%centers_set(ispin)%array,2)
653 : ! IF (nmoloc/=nmo) CALL cp_abort(__LOCATION__,&
654 : ! "The nbr of localized functions "//&
655 : ! & "is not equal to the nbr of states.")
656 : !
657 : ! check that the occupation is uniform
658 674 : IF (.NOT. uniform_occupation) CPABORT("nmo != homo")
659 : END DO
660 : !
661 : !
662 174 : logger => cp_get_default_logger()
663 174 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
664 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
665 174 : extension=".linresLog")
666 :
667 174 : IF (output_unit > 0) THEN
668 87 : WRITE (output_unit, "(/,T20,A,/)") "*** Start current Calculation ***"
669 87 : WRITE (output_unit, "(T10,A,/)") "Inizialization of the current environment"
670 : END IF
671 :
672 174 : CALL current_env_cleanup(current_env)
673 174 : current_env%gauge_init = .FALSE.
674 4698 : current_env%chi_tensor(:, :, :) = 0.0_dp
675 4698 : current_env%chi_tensor_loc(:, :, :) = 0.0_dp
676 174 : current_env%nao = nao
677 174 : current_env%full = .TRUE.
678 174 : current_env%do_selected_states = .FALSE.
679 : !
680 : ! If current_density or full_nmr different allocations are required
681 174 : current_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%CURRENT")
682 : !
683 : ! Select the gauge
684 174 : CALL section_vals_val_get(current_section, "GAUGE", i_val=current_env%gauge)
685 178 : SELECT CASE (current_env%gauge)
686 : CASE (current_gauge_r)
687 4 : current_env%gauge_name = "R"
688 : CASE (current_gauge_r_and_step_func)
689 142 : current_env%gauge_name = "R_AND_STEP_FUNCTION"
690 : CASE (current_gauge_atom)
691 28 : current_env%gauge_name = "ATOM"
692 : CASE DEFAULT
693 174 : CPABORT("Unknown gauge, try again...")
694 : END SELECT
695 : !
696 : ! maximal radius to build the atom gauge
697 : CALL section_vals_val_get(current_section, "GAUGE_ATOM_RADIUS", &
698 174 : r_val=current_env%gauge_atom_radius)
699 : ! use old gauge atom
700 : CALL section_vals_val_get(current_section, "USE_OLD_GAUGE_ATOM", &
701 174 : l_val=current_env%use_old_gauge_atom)
702 : ! chi for pbc
703 174 : CALL section_vals_val_get(current_section, "CHI_PBC", l_val=current_env%chi_pbc)
704 : !
705 : ! use old gauge atom
706 : CALL section_vals_val_get(current_section, "USE_OLD_GAUGE_ATOM", &
707 174 : l_val=current_env%use_old_gauge_atom)
708 : !
709 : ! chi for pbc
710 174 : CALL section_vals_val_get(current_section, "CHI_PBC", l_val=current_env%chi_pbc)
711 : !
712 : ! which center for the orbitals shall we use
713 174 : CALL section_vals_val_get(current_section, "ORBITAL_CENTER", i_val=current_env%orb_center)
714 300 : SELECT CASE (current_env%orb_center)
715 : CASE (current_orb_center_wannier)
716 : !
717 126 : current_env%orb_center_name = "WANNIER"
718 : CASE (current_orb_center_common)
719 : !
720 32 : current_env%orb_center_name = "COMMON"
721 32 : current_env%full = .FALSE.
722 : !
723 : ! Is there a user specified common_center?
724 32 : CALL section_vals_val_get(current_section, "COMMON_CENTER", r_vals=common_center)
725 : CASE (current_orb_center_atom)
726 : !
727 10 : current_env%orb_center_name = "ATOM"
728 : CASE (current_orb_center_box)
729 : !
730 6 : current_env%orb_center_name = "BOX"
731 : !
732 : ! Is there a user specified nbox?
733 6 : CALL section_vals_val_get(current_section, "NBOX", i_vals=nbox)
734 : CASE DEFAULT
735 174 : CPABORT("Unknown orbital center, try again...")
736 : END SELECT
737 :
738 174 : CALL section_vals_val_get(current_section, "FORCE_NO_FULL", l_val=force_no_full)
739 174 : IF (force_no_full) current_env%full = .FALSE.
740 : !
741 : ! Check if restat also psi0 should be restarted
742 : !IF(current_env%restart_current .AND. scf_control%density_guess/=restart_guess)THEN
743 : ! CPABORT("restart_nmr requires density_guess=restart")
744 : !ENDIF
745 : !
746 : ! check that the psi0 are localized and you have all the centers
747 174 : CPASSERT(linres_control%localized_psi0)
748 174 : IF (output_unit > 0) THEN
749 : WRITE (output_unit, '(A)') &
750 87 : ' To get CURRENT parameters within PBC you need localized zero order orbitals '
751 : END IF
752 174 : qs_loc_env => linres_control%qs_loc_env
753 174 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
754 : !
755 : !
756 : ALLOCATE (current_env%centers_set(nspins), current_env%center_list(nspins), &
757 1718 : state_list(max_states, nspins))
758 2196 : state_list(:, :) = HUGE(0)
759 522 : nstate_list(:) = HUGE(0)
760 : !
761 : !
762 : !
763 : ! if qmmm is selected change the definition of the center of the box 0 -> L/2
764 174 : IF (current_env%do_qmmm) THEN
765 12 : DO ispin = 1, nspins
766 42 : DO istate = 1, SIZE(localized_wfn_control%centers_set(ispin)%array, 2)
767 : ! just to be sure...
768 30 : r(:) = pbc(localized_wfn_control%centers_set(ispin)%array(:, istate), cell)
769 30 : IF (r(1) < 0.0_dp) THEN
770 : localized_wfn_control%centers_set(ispin)%array(1, istate) = &
771 6 : r(1) + cell%hmat(1, 1)
772 : END IF
773 30 : IF (r(2) < 0.0_dp) THEN
774 : localized_wfn_control%centers_set(ispin)%array(2, istate) = &
775 24 : r(2) + cell%hmat(2, 2)
776 : END IF
777 36 : IF (r(3) < 0.0_dp) THEN
778 : localized_wfn_control%centers_set(ispin)%array(3, istate) = &
779 6 : r(3) + cell%hmat(3, 3)
780 : END IF
781 : END DO
782 : END DO
783 : END IF
784 : !
785 : !
786 : !
787 : ! if the user has requested to compute the response for a subset of the states
788 : ! we collect them here. it requies the states to be localized!
789 : CALL section_vals_val_get(current_section, "SELECTED_STATES_ATOM_RADIUS", &
790 174 : r_val=current_env%selected_states_atom_radius)
791 174 : CALL section_vals_val_get(current_section, "SELECTED_STATES_ON_ATOM_LIST", n_rep_val=n_rep)
792 : !
793 174 : current_env%do_selected_states = n_rep > 0
794 : !
795 : ! for the moment selected states doesnt work with the preconditioner FULL_ALL
796 174 : IF (linres_control%preconditioner_type == ot_precond_full_all .AND. current_env%do_selected_states) THEN
797 0 : CPABORT("Selected states doesnt work with the preconditioner FULL_ALL")
798 : END IF
799 : !
800 : !
801 174 : NULLIFY (current_env%selected_states_on_atom_list)
802 174 : n = 0
803 182 : DO ir = 1, n_rep
804 8 : NULLIFY (list)
805 : CALL section_vals_val_get(current_section, "SELECTED_STATES_ON_ATOM_LIST", &
806 8 : i_rep_val=ir, i_vals=list)
807 182 : IF (ASSOCIATED(list)) THEN
808 8 : CALL reallocate(current_env%selected_states_on_atom_list, 1, n + SIZE(list))
809 32 : DO ini = 1, SIZE(list)
810 32 : current_env%selected_states_on_atom_list(ini + n) = list(ini)
811 : END DO
812 8 : n = n + SIZE(list)
813 : END IF
814 : END DO
815 : !
816 : ! build the subset
817 174 : IF (current_env%do_selected_states) THEN
818 8 : selected_states_on_atom_list => current_env%selected_states_on_atom_list
819 20 : DO ispin = 1, nspins
820 12 : center_array => localized_wfn_control%centers_set(ispin)%array
821 12 : nstate = 0
822 184 : DO istate = 1, SIZE(center_array, 2)
823 436 : DO i = 1, SIZE(selected_states_on_atom_list, 1)
824 364 : iatom = selected_states_on_atom_list(i)
825 1456 : r(:) = pbc(center_array(1:3, istate) - particle_set(iatom)%r(:), cell)
826 : ! SQRT(DOT_PRODUCT(r, r)) <= current_env%selected_states_atom_radius
827 1456 : IF ((DOT_PRODUCT(r, r)) <= (current_env%selected_states_atom_radius &
828 : *current_env%selected_states_atom_radius)) &
829 60 : THEN
830 : !
831 : ! add the state to the list
832 112 : nstate = nstate + 1
833 112 : state_list(nstate, ispin) = istate
834 112 : EXIT
835 : END IF
836 : END DO
837 : END DO
838 20 : nstate_list(ispin) = nstate
839 : END DO
840 : ELSE
841 404 : DO ispin = 1, nspins
842 238 : center_array => localized_wfn_control%centers_set(ispin)%array
843 238 : nstate = 0
844 1724 : DO istate = 1, SIZE(center_array, 2)
845 1486 : nstate = nstate + 1
846 1724 : state_list(nstate, ispin) = istate
847 : END DO
848 404 : nstate_list(ispin) = nstate
849 : END DO
850 : END IF
851 : !
852 : !
853 : !
854 : ! clustering the states
855 424 : DO ispin = 1, nspins
856 250 : nstate = nstate_list(ispin)
857 250 : current_env%nstates(ispin) = nstate
858 : !
859 : ALLOCATE (current_env%center_list(ispin)%array(2, nstate + 1), &
860 1250 : current_env%centers_set(ispin)%array(3, nstate))
861 5794 : current_env%center_list(ispin)%array(:, :) = HUGE(0)
862 6642 : current_env%centers_set(ispin)%array(:, :) = HUGE(0.0_dp)
863 : !
864 250 : center_array => localized_wfn_control%centers_set(ispin)%array
865 : !
866 : ! point to the psi0 centers
867 174 : SELECT CASE (current_env%orb_center)
868 : CASE (current_orb_center_wannier)
869 : !
870 : ! use the wannier center as -center-
871 180 : current_env%nbr_center(ispin) = nstate
872 1284 : DO is = 1, nstate
873 1104 : istate = state_list(is, ispin)
874 4416 : current_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
875 1104 : current_env%center_list(ispin)%array(1, is) = is
876 1284 : current_env%center_list(ispin)%array(2, is) = istate
877 : END DO
878 180 : current_env%center_list(ispin)%array(1, nstate + 1) = nstate + 1
879 : !
880 : CASE (current_orb_center_common)
881 : !
882 : ! use a common -center-
883 176 : current_env%centers_set(ispin)%array(:, 1) = common_center(:)
884 44 : current_env%nbr_center(ispin) = 1
885 44 : current_env%center_list(ispin)%array(1, 1) = 1
886 44 : current_env%center_list(ispin)%array(1, 2) = nstate + 1
887 310 : DO is = 1, nstate
888 266 : istate = state_list(is, ispin)
889 310 : current_env%center_list(ispin)%array(2, is) = istate
890 : END DO
891 : !
892 : CASE (current_orb_center_atom)
893 : !
894 : ! use the atom as -center-
895 48 : ALLOCATE (buff(nstate_list(ispin)))
896 16 : buff(:) = 0
897 : !
898 158 : DO is = 1, nstate
899 142 : istate = state_list(is, ispin)
900 142 : mdist = HUGE(0.0_dp)
901 712 : DO iatom = 1, natom
902 554 : r = pbc(particle_set(iatom)%r(:), cell)
903 554 : rab = pbc(r, center_array(1:3, istate), cell)
904 554 : dist = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
905 696 : IF (dist < mdist) THEN
906 322 : buff(is) = iatom
907 322 : mdist = dist
908 : END IF
909 : END DO
910 : END DO
911 : !
912 16 : i = 0
913 16 : ii = 1
914 16 : current_env%center_list(ispin)%array(1, 1) = 1
915 76 : DO iatom = 1, natom
916 60 : j = 0
917 60 : is0 = .TRUE.
918 614 : DO is = 1, nstate
919 554 : istate = state_list(is, ispin)
920 614 : IF (buff(is) == iatom) THEN
921 142 : j = j + 1
922 142 : i = i + 1
923 142 : is0 = .FALSE.
924 142 : current_env%center_list(ispin)%array(2, i) = istate
925 : END IF
926 : END DO
927 76 : IF (.NOT. is0) THEN
928 48 : IF (output_unit > 0) THEN
929 24 : WRITE (output_unit, '(T2,A,I6,A,I6)') 'clustering ', j, ' center(s) on atom ', iatom
930 : END IF
931 : current_env%center_list(ispin)%array(1, ii + 1) = &
932 48 : current_env%center_list(ispin)%array(1, ii) + j
933 : current_env%centers_set(ispin)%array(:, ii) = &
934 192 : pbc(particle_set(iatom)%r, cell)
935 48 : ii = ii + 1
936 : END IF
937 : END DO
938 16 : current_env%nbr_center(ispin) = ii - 1
939 : !
940 16 : DEALLOCATE (buff)
941 : CASE (current_orb_center_box)
942 : !
943 : ! use boxes as -center-
944 10 : nbr_box = nbox(1)*nbox(2)*nbox(3)
945 50 : ALLOCATE (rbuff(3, nbr_box), buff(nstate))
946 3546 : rbuff(:, :) = HUGE(0.0_dp)
947 10 : buff(:) = 0
948 : !
949 10 : ibox = 1
950 54 : DO iz = 1, nbox(3)
951 250 : DO iy = 1, nbox(2)
952 1124 : DO ix = 1, nbox(1)
953 884 : rbuff(1, ibox) = cell%hmat(1, 1)*((REAL(ix, dp) - 0.5_dp)/REAL(nbox(1), dp) - 0.5_dp)
954 884 : rbuff(2, ibox) = cell%hmat(2, 2)*((REAL(iy, dp) - 0.5_dp)/REAL(nbox(2), dp) - 0.5_dp)
955 884 : rbuff(3, ibox) = cell%hmat(3, 3)*((REAL(iz, dp) - 0.5_dp)/REAL(nbox(3), dp) - 0.5_dp)
956 1080 : ibox = ibox + 1
957 : END DO
958 : END DO
959 : END DO
960 : !
961 96 : DO is = 1, nstate
962 86 : istate = state_list(is, ispin)
963 86 : mdist = HUGE(0.0_dp)
964 8406 : DO ibox = 1, nbr_box
965 8310 : rab(:) = pbc(rbuff(:, ibox), center_array(1:3, istate), cell)
966 8310 : dist = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
967 8396 : IF (dist < mdist) THEN
968 628 : buff(is) = ibox
969 628 : mdist = dist
970 : END IF
971 : END DO
972 : END DO
973 : !
974 10 : i = 0
975 10 : ii = 1
976 10 : current_env%center_list(ispin)%array(1, 1) = 1
977 894 : DO ibox = 1, nbr_box
978 884 : j = 0
979 884 : is0 = .TRUE.
980 9194 : DO is = 1, nstate
981 8310 : istate = state_list(is, ispin)
982 9194 : IF (buff(is) == ibox) THEN
983 86 : j = j + 1
984 86 : i = i + 1
985 86 : is0 = .FALSE.
986 86 : current_env%center_list(ispin)%array(2, i) = istate
987 : END IF
988 : END DO
989 894 : IF (.NOT. is0) THEN
990 54 : IF (output_unit > 0) THEN
991 27 : WRITE (output_unit, '(T2,A,I6,A,I6)') 'clustering ', j, ' center(s) on box ', ibox
992 : END IF
993 : current_env%center_list(ispin)%array(1, ii + 1) = &
994 54 : current_env%center_list(ispin)%array(1, ii) + j
995 216 : current_env%centers_set(ispin)%array(:, ii) = rbuff(:, ibox)
996 : ii = ii + 1
997 : END IF
998 : END DO
999 10 : current_env%nbr_center(ispin) = ii - 1
1000 : !
1001 10 : DEALLOCATE (buff, rbuff)
1002 : CASE DEFAULT
1003 250 : CPABORT("Unknown orbital center, try again...")
1004 : END SELECT
1005 : END DO
1006 : !
1007 : !
1008 : !
1009 : ! block the states for each center
1010 772 : ALLOCATE (current_env%psi0_order(nspins))
1011 424 : DO ispin = 1, nspins
1012 250 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1013 250 : NULLIFY (tmp_fm_struct)
1014 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
1015 : ncol_global=nstate_list(ispin), para_env=para_env, &
1016 250 : context=mo_coeff%matrix_struct%context)
1017 250 : CALL cp_fm_create(current_env%psi0_order(ispin), tmp_fm_struct)
1018 250 : CALL cp_fm_struct_release(tmp_fm_struct)
1019 424 : CALL cp_fm_set_all(current_env%psi0_order(ispin), 0.0_dp)
1020 : END DO
1021 : !
1022 : !
1023 424 : DO ispin = 1, nspins
1024 250 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1025 250 : jstate = 0
1026 1674 : DO icenter = 1, current_env%nbr_center(ispin)
1027 2848 : DO i = current_env%center_list(ispin)%array(1, icenter), &
1028 1500 : current_env%center_list(ispin)%array(1, icenter + 1) - 1
1029 : !
1030 1598 : istate = current_env%center_list(ispin)%array(2, i)
1031 1598 : jstate = jstate + 1
1032 : !
1033 2848 : IF (current_env%do_selected_states) THEN ! this should be removed. always reorder the states
1034 : ! the blocking works (so far) with all the precond except FULL_ALL
1035 : ! the center_list(ispin)%array(2,i) array is obsolete then
1036 112 : current_env%center_list(ispin)%array(2, i) = jstate
1037 112 : CALL cp_fm_to_fm(mo_coeff, current_env%psi0_order(ispin), 1, istate, jstate)
1038 : ELSE
1039 : ! for the moment we just copy them
1040 1486 : CALL cp_fm_to_fm(mo_coeff, current_env%psi0_order(ispin), 1, istate, istate)
1041 : END IF
1042 : END DO
1043 : END DO
1044 : END DO
1045 : !
1046 : !
1047 : !
1048 : !
1049 174 : IF (current_env%do_qmmm .AND.&
1050 : & (current_env%orb_center == current_orb_center_wannier .OR. &
1051 : current_env%orb_center == current_orb_center_atom .OR. &
1052 : current_env%orb_center == current_orb_center_box)) THEN
1053 4 : IF (output_unit > 0) THEN
1054 2 : WRITE (output_unit, *) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING '
1055 2 : WRITE (output_unit, *) 'orbital center shifted to match the '
1056 2 : WRITE (output_unit, *) 'center of the box (L/2 L/2 L/2) (superdirty...)'
1057 2 : WRITE (output_unit, *) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING '
1058 : END IF
1059 8 : DO ispin = 1, nspins
1060 24 : DO istate = 1, current_env%nbr_center(ispin)
1061 16 : IF (current_env%centers_set(ispin)%array(1, istate) <= 0.0_dp) THEN
1062 : current_env%centers_set(ispin)%array(1, istate) = &
1063 0 : current_env%centers_set(ispin)%array(1, istate) + cell%hmat(1, 1)
1064 : END IF
1065 16 : IF (current_env%centers_set(ispin)%array(2, istate) <= 0.0_dp) THEN
1066 : current_env%centers_set(ispin)%array(2, istate) = &
1067 0 : current_env%centers_set(ispin)%array(2, istate) + cell%hmat(2, 2)
1068 : END IF
1069 20 : IF (current_env%centers_set(ispin)%array(3, istate) <= 0.0_dp) THEN
1070 : current_env%centers_set(ispin)%array(3, istate) = &
1071 0 : current_env%centers_set(ispin)%array(3, istate) + cell%hmat(3, 3)
1072 : END IF
1073 : END DO
1074 : END DO
1075 : ! printing
1076 8 : DO ispin = 1, nspins
1077 178 : IF (output_unit > 0) THEN
1078 2 : WRITE (output_unit, '(/,T2,A,I2)') "WANNIER CENTERS for spin ", ispin
1079 2 : WRITE (output_unit, '(/,T18,A)') "--------------- Shifted Centers --------------- "
1080 10 : DO istate = 1, current_env%nbr_center(ispin)
1081 : WRITE (output_unit, '(T5,A,I6,3F12.6)') &
1082 34 : 'state ', istate, current_env%centers_set(ispin)%array(1:3, istate)
1083 : END DO
1084 : END IF
1085 : END DO
1086 : END IF
1087 : !
1088 : !
1089 : !
1090 : ! reorder the center dependent responses
1091 424 : max_nbr_center = MAXVAL(current_env%nbr_center(1:nspins))
1092 174 : current_env%nao = nao
1093 696 : ALLOCATE (current_env%statetrueindex(3, max_nbr_center, nspins))
1094 5792 : current_env%statetrueindex(:, :, :) = 0
1095 : IF (.TRUE.) THEN
1096 522 : ALLOCATE (state_done(3, max_nbr_center))
1097 424 : DO ispin = 1, nspins
1098 250 : state_done(:, :) = .FALSE.
1099 250 : current_env%statetrueindex(1, 1, ispin) = 1
1100 250 : center(1) = current_env%centers_set(ispin)%array(1, 1)
1101 250 : center(2) = current_env%centers_set(ispin)%array(2, 1)
1102 250 : center(3) = current_env%centers_set(ispin)%array(3, 1)
1103 250 : state_done(1, 1) = .TRUE.
1104 250 : icount = 1
1105 250 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
1106 : !
1107 1174 : DO idir = 1, 3
1108 750 : ini = 1
1109 750 : IF (idir == 1) ini = 2
1110 4250 : DO istate = ini, current_env%nbr_center(ispin)
1111 3500 : mdist = HUGE(0.0_dp)
1112 : !
1113 26496 : DO istate2 = 1, current_env%nbr_center(ispin)
1114 26496 : IF (.NOT. state_done(idir, istate2)) THEN
1115 12748 : center2(1) = current_env%centers_set(ispin)%array(1, istate2)
1116 12748 : center2(2) = current_env%centers_set(ispin)%array(2, istate2)
1117 12748 : center2(3) = current_env%centers_set(ispin)%array(3, istate2)
1118 : !
1119 12748 : rab = pbc(center, center2, cell)
1120 12748 : CALL set_vecp(idir, j, k)
1121 12748 : dist = SQRT(rab(j)*rab(j) + rab(k)*rab(k))
1122 : !
1123 12748 : IF (dist < mdist) THEN
1124 6644 : mdist = dist
1125 6644 : istate_next = istate2
1126 : END IF
1127 : END IF
1128 : END DO ! istate2
1129 : !
1130 3500 : icount = icount + 1
1131 3500 : state_done(idir, istate_next) = .TRUE.
1132 3500 : current_env%statetrueindex(idir, icount, ispin) = istate_next
1133 3500 : center(1) = current_env%centers_set(ispin)%array(1, istate_next)
1134 3500 : center(2) = current_env%centers_set(ispin)%array(2, istate_next)
1135 4250 : center(3) = current_env%centers_set(ispin)%array(3, istate_next)
1136 : END DO ! istate
1137 1000 : icount = 0
1138 : END DO ! idir
1139 : END DO
1140 174 : DEALLOCATE (state_done)
1141 : ELSE
1142 : DO ispin = 1, nspins
1143 : DO icenter = 1, current_env%nbr_center(ispin)
1144 : current_env%statetrueindex(:, icenter, ispin) = icenter
1145 : END DO
1146 : END DO
1147 : END IF
1148 : !
1149 : !
1150 174 : IF (output_unit > 0) THEN
1151 87 : WRITE (output_unit, "(T2,A,T60,A)") "CURRENT| Gauge used", &
1152 507 : REPEAT(' ', 20 - LEN_TRIM(current_env%gauge_name))//TRIM(current_env%gauge_name)
1153 87 : WRITE (output_unit, "(T2,A,T79,L1)") "CURRENT| Use old gauge code", current_env%use_old_gauge_atom
1154 87 : WRITE (output_unit, "(T2,A,T79,L1)") "CURRENT| Compute chi for PBC calculation", current_env%chi_pbc
1155 87 : IF (current_env%gauge == current_gauge_atom) THEN
1156 14 : WRITE (output_unit, "(T2,A,T68,E12.6)") "CURRENT| Radius for building the gauge", &
1157 28 : current_env%gauge_atom_radius
1158 : END IF
1159 87 : WRITE (output_unit, "(T2,A,T60,A)") "CURRENT| Orbital center used", &
1160 1348 : REPEAT(' ', 20 - LEN_TRIM(current_env%orb_center_name))//TRIM(current_env%orb_center_name)
1161 87 : IF (current_env%orb_center == current_orb_center_common) THEN
1162 64 : WRITE (output_unit, "(T2,A,T50,3F10.6)") "CURRENT| Common center", common_center(1:3)
1163 71 : ELSE IF (current_env%orb_center == current_orb_center_box) THEN
1164 12 : WRITE (output_unit, "(T2,A,T60,3I5)") "CURRENT| Nbr boxes in each direction", nbox(1:3)
1165 : END IF
1166 : !
1167 212 : DO ispin = 1, nspins
1168 125 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
1169 125 : WRITE (output_unit, '(T2,A,I6,A,I6,A,I2)') 'CURRENT| Compute', &
1170 125 : nstate_list(ispin), ' selected response functions out of', &
1171 250 : nmo, ' for spin ', ispin
1172 : !
1173 125 : WRITE (output_unit, '(T2,A,I6,A,I2)') 'CURRENT| There is a total of ', &
1174 462 : current_env%nbr_center(ispin), ' (clustered) center(s) for spin ', ispin
1175 : END DO
1176 : END IF
1177 :
1178 174 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
1179 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
1180 :
1181 0 : NULLIFY (bounds, list)
1182 0 : ncubes = 0
1183 : CALL section_vals_val_get(current_section,&
1184 : & "PRINT%RESPONSE_FUNCTION_CUBES%CUBES_LU_BOUNDS", &
1185 0 : i_vals=bounds)
1186 0 : ncubes = bounds(2) - bounds(1) + 1
1187 0 : IF (ncubes > 0) THEN
1188 0 : ALLOCATE (current_env%list_cubes(ncubes))
1189 0 : DO ir = 1, ncubes
1190 0 : current_env%list_cubes(ir) = bounds(1) + (ir - 1)
1191 : END DO
1192 : END IF
1193 0 : IF (.NOT. ASSOCIATED(current_env%list_cubes)) THEN
1194 : CALL section_vals_val_get(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%CUBES_LIST", &
1195 0 : n_rep_val=n_rep)
1196 0 : ncubes = 0
1197 0 : DO ir = 1, n_rep
1198 0 : NULLIFY (list)
1199 : CALL section_vals_val_get(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%CUBES_LIST", &
1200 0 : i_rep_val=ir, i_vals=list)
1201 0 : IF (ASSOCIATED(list)) THEN
1202 0 : CALL reallocate(current_env%list_cubes, 1, ncubes + SIZE(list))
1203 0 : DO ini = 1, SIZE(list)
1204 0 : current_env%list_cubes(ini + ncubes) = list(ini)
1205 : END DO
1206 0 : ncubes = ncubes + SIZE(list)
1207 : END IF
1208 : END DO ! ir
1209 : END IF
1210 0 : IF (.NOT. ASSOCIATED(current_env%list_cubes)) THEN
1211 0 : ALLOCATE (current_env%list_cubes(max_states))
1212 0 : DO ir = 1, max_states
1213 0 : current_env%list_cubes(ir) = ir
1214 : END DO
1215 : END IF
1216 : END IF
1217 174 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
1218 : & "PRINT%CURRENT_CUBES"), cp_p_file)) THEN
1219 : END IF
1220 :
1221 : ! for the chemical shift we need 6 psi1, i.e. 6 optimization procedures
1222 : ! They become 9 if full nmr is calculated, i.e. with the correction term too
1223 : ! All of them are required at the end of the optimization procedure
1224 : ! if the current density and the induced field have to be calculated
1225 : ! If instead only the shift is needed, only one psi1 should be enough, providing
1226 : ! that after every optimization the corresponding shift contribution is calculated
1227 : ! prepare the psi1
1228 : !
1229 : ! for the rxp we cannot calculate it a priori because it is in facts (r-dk)xp
1230 : ! where dk is the center of the orbital it is applied to. We would need nstate operators
1231 : ! What we can store is (r-dk)xp|psi0k> for each k, which is a full matrix only
1232 : ! Therefore we prepare here the full matrix p_psi0 and rxp_psi0
1233 : ! We also need a temporary sparse matrix where to store the integrals during the calculation
1234 : ALLOCATE (current_env%p_psi0(nspins, 3), current_env%rxp_psi0(nspins, 3), &
1235 6132 : current_env%psi1_p(nspins, 3), current_env%psi1_rxp(nspins, 3))
1236 174 : IF (current_env%full) THEN
1237 1328 : ALLOCATE (current_env%psi1_D(nspins, 3))
1238 : END IF
1239 424 : DO ispin = 1, nspins
1240 250 : mo_coeff => current_env%psi0_order(ispin)
1241 250 : NULLIFY (tmp_fm_struct)
1242 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
1243 : ncol_global=current_env%nstates(ispin), &
1244 250 : context=mo_coeff%matrix_struct%context)
1245 1000 : DO idir = 1, 3
1246 750 : CALL cp_fm_create(current_env%psi1_p(ispin, idir), tmp_fm_struct)
1247 750 : CALL cp_fm_create(current_env%psi1_rxp(ispin, idir), tmp_fm_struct)
1248 750 : CALL cp_fm_create(current_env%p_psi0(ispin, idir), tmp_fm_struct)
1249 750 : CALL cp_fm_create(current_env%rxp_psi0(ispin, idir), tmp_fm_struct)
1250 1000 : IF (current_env%full) THEN
1251 618 : CALL cp_fm_create(current_env%psi1_D(ispin, idir), tmp_fm_struct)
1252 : END IF
1253 : END DO
1254 424 : CALL cp_fm_struct_release(tmp_fm_struct)
1255 : END DO
1256 : !
1257 : ! If the current density on the grid needs to be stored
1258 174 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1259 696 : ALLOCATE (current_env%jrho1_set(3))
1260 696 : DO idir = 1, 3
1261 522 : NULLIFY (rho_r, rho_g)
1262 4110 : ALLOCATE (rho_r(nspins), rho_g(nspins))
1263 1272 : DO ispin = 1, nspins
1264 750 : CALL auxbas_pw_pool%create_pw(rho_r(ispin))
1265 750 : CALL pw_zero(rho_r(ispin))
1266 750 : CALL auxbas_pw_pool%create_pw(rho_g(ispin))
1267 1272 : CALL pw_zero(rho_g(ispin))
1268 : END DO
1269 522 : NULLIFY (current_env%jrho1_set(idir)%rho)
1270 522 : ALLOCATE (current_env%jrho1_set(idir)%rho)
1271 522 : CALL qs_rho_create(current_env%jrho1_set(idir)%rho)
1272 : CALL qs_rho_set(current_env%jrho1_set(idir)%rho, &
1273 696 : rho_r=rho_r, rho_g=rho_g)
1274 : END DO
1275 : !
1276 : ! Initialize local current density if GAPW calculation
1277 174 : IF (gapw) THEN
1278 96 : CALL init_jrho_atom_set(jrho1_atom_set, atomic_kind_set, nspins)
1279 96 : CALL set_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set)
1280 : END IF
1281 : !
1282 : ALLOCATE (current_env%basisfun_center(3, current_env%nao), &
1283 1044 : first_sgf(natom), last_sgf(natom))
1284 15886 : current_env%basisfun_center = 0.0_dp
1285 : CALL get_particle_set(particle_set, qs_kind_set, &
1286 : first_sgf=first_sgf, &
1287 174 : last_sgf=last_sgf)
1288 : !Build 3 arrays where for each contracted basis function
1289 : !the x y and z coordinates of the center are given
1290 710 : DO iatom = 1, natom
1291 4638 : DO iao = first_sgf(iatom), last_sgf(iatom)
1292 16248 : DO idir = 1, 3
1293 15712 : current_env%basisfun_center(idir, iao) = particle_set(iatom)%r(idir)
1294 : END DO
1295 : END DO
1296 : END DO
1297 :
1298 174 : DEALLOCATE (first_sgf, last_sgf, state_list)
1299 :
1300 1218 : current_env%simple_done(1:6) = .FALSE.
1301 :
1302 696 : ALLOCATE (current_env%full_done(3*nspins, max_nbr_center))
1303 5052 : current_env%full_done = .FALSE.
1304 : !
1305 : ! allocate pointer for the gauge
1306 5900 : ALLOCATE (current_env%rs_gauge(3, pw_env%gridlevel_info%ngrid_levels))
1307 :
1308 174 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, "PRINT%PROGRAM_RUN_INFO")
1309 174 : CALL timestop(handle)
1310 :
1311 870 : END SUBROUTINE current_env_init
1312 :
1313 : ! **************************************************************************************************
1314 : !> \brief ...
1315 : !> \param current_env ...
1316 : ! **************************************************************************************************
1317 348 : SUBROUTINE current_env_cleanup(current_env)
1318 :
1319 : TYPE(current_env_type) :: current_env
1320 :
1321 : INTEGER :: i, idir, ispin
1322 :
1323 348 : CALL cp_fm_release(current_env%psi0_order)
1324 : !
1325 : !psi1_p
1326 348 : CALL cp_fm_release(current_env%psi1_p)
1327 : !
1328 : !psi1_rxp
1329 348 : CALL cp_fm_release(current_env%psi1_rxp)
1330 : !
1331 : !psi1_D
1332 348 : CALL cp_fm_release(current_env%psi1_D)
1333 : !
1334 : !p_psi0
1335 348 : CALL cp_fm_release(current_env%p_psi0)
1336 : !
1337 : !rxp_psi0
1338 348 : CALL cp_fm_release(current_env%rxp_psi0)
1339 : !
1340 348 : IF (ASSOCIATED(current_env%centers_set)) THEN
1341 424 : DO ispin = 1, SIZE(current_env%centers_set, 1)
1342 424 : DEALLOCATE (current_env%centers_set(ispin)%array)
1343 : END DO
1344 174 : DEALLOCATE (current_env%centers_set)
1345 : END IF
1346 :
1347 348 : IF (ASSOCIATED(current_env%center_list)) THEN
1348 424 : DO ispin = 1, SIZE(current_env%center_list, 1)
1349 424 : DEALLOCATE (current_env%center_list(ispin)%array)
1350 : END DO
1351 174 : DEALLOCATE (current_env%center_list)
1352 : END IF
1353 :
1354 348 : IF (ASSOCIATED(current_env%list_cubes)) THEN
1355 0 : DEALLOCATE (current_env%list_cubes)
1356 : END IF
1357 : ! Current density on the grid
1358 348 : IF (ASSOCIATED(current_env%jrho1_set)) THEN
1359 696 : DO idir = 1, 3
1360 522 : CALL qs_rho_clear(current_env%jrho1_set(idir)%rho)
1361 696 : DEALLOCATE (current_env%jrho1_set(idir)%rho)
1362 : END DO
1363 174 : DEALLOCATE (current_env%jrho1_set)
1364 : END IF
1365 : !
1366 : ! Local current density, atom by atom (only gapw)
1367 348 : IF (ASSOCIATED(current_env%jrho1_atom_set)) THEN
1368 96 : CALL deallocate_jrho_atom_set(current_env%jrho1_atom_set)
1369 : END IF
1370 :
1371 : !fullnmr_done
1372 348 : IF (ASSOCIATED(current_env%full_done)) THEN
1373 174 : DEALLOCATE (current_env%full_done)
1374 : END IF
1375 348 : IF (ASSOCIATED(current_env%basisfun_center)) THEN
1376 174 : DEALLOCATE (current_env%basisfun_center)
1377 : END IF
1378 348 : IF (ASSOCIATED(current_env%statetrueindex)) THEN
1379 174 : DEALLOCATE (current_env%statetrueindex)
1380 : END IF
1381 348 : IF (ASSOCIATED(current_env%selected_states_on_atom_list)) THEN
1382 8 : DEALLOCATE (current_env%selected_states_on_atom_list)
1383 : END IF
1384 : ! give back the gauge
1385 348 : IF (ASSOCIATED(current_env%rs_gauge)) THEN
1386 696 : DO idir = 1, 3
1387 2772 : DO i = 1, SIZE(current_env%rs_gauge, 2)
1388 2598 : CALL rs_grid_release(current_env%rs_gauge(idir, i))
1389 : END DO
1390 : END DO
1391 174 : DEALLOCATE (current_env%rs_gauge)
1392 : END IF
1393 :
1394 : ! give back the buf
1395 348 : IF (ASSOCIATED(current_env%rs_buf)) THEN
1396 136 : DO i = 1, SIZE(current_env%rs_buf)
1397 136 : CALL rs_grid_release(current_env%rs_buf(i))
1398 : END DO
1399 28 : DEALLOCATE (current_env%rs_buf)
1400 : END IF
1401 :
1402 348 : END SUBROUTINE current_env_cleanup
1403 :
1404 : END MODULE qs_linres_current_utils
|
