Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
30 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
31 : dbcsr_allocate_matrix_set,&
32 : dbcsr_deallocate_matrix_set
33 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
34 : cp_fm_scale_and_add
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_release,&
37 : cp_fm_struct_type
38 : USE cp_fm_types, ONLY: cp_fm_create,&
39 : cp_fm_release,&
40 : cp_fm_set_all,&
41 : cp_fm_to_fm,&
42 : cp_fm_type
43 : USE cp_log_handling, ONLY: cp_get_default_logger,&
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_should_output,&
49 : cp_print_key_unit_nr
50 : USE dbcsr_api, ONLY: dbcsr_convert_offsets_to_sizes,&
51 : dbcsr_copy,&
52 : dbcsr_create,&
53 : dbcsr_distribution_type,&
54 : dbcsr_p_type,&
55 : dbcsr_set,&
56 : dbcsr_type_antisymmetric
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_path_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_path_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 : nze=0, 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 : !CALL p_xyz_ao(op_p_ao,qs_env,minimum_image=.FALSE.)
263 142 : CALL build_lin_mom_matrix(qs_env, op_p_ao)
264 : !
265 : END IF
266 : !
267 : ! set response to zero
268 : !
269 696 : DO idir = 1, 3
270 1446 : DO ispin = 1, nspins
271 750 : CALL cp_fm_set_all(psi1_p(ispin, idir), 0.0_dp)
272 750 : CALL cp_fm_set_all(psi1_rxp(ispin, idir), 0.0_dp)
273 1272 : IF (current_env%full) CALL cp_fm_set_all(psi1_D(ispin, idir), 0.0_dp)
274 : END DO
275 : END DO
276 : !
277 : !
278 : !
279 : ! load restart file
280 : !
281 174 : first_center = 0
282 174 : IF (linres_control%linres_restart) THEN
283 24 : DO idir = 1, 3
284 : ! operator p
285 18 : CALL linres_read_restart(qs_env, lr_section, psi1_p(:, idir), idir, "nmr_p")
286 : ! operator rxp
287 24 : CALL linres_read_restart(qs_env, lr_section, psi1_rxp(:, idir), idir, "nmr_rxp")
288 : END DO
289 : END IF
290 : !
291 : !
292 : !
293 : !
294 174 : should_stop = .FALSE.
295 174 : current_env%all_pert_op_done = .FALSE.
296 : ! operator p
297 696 : DO idir = 1, 3
298 522 : IF (should_stop) EXIT
299 522 : IF (output_unit > 0) THEN
300 261 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator P_"//ACHAR(idir + 119)
301 : END IF
302 : !
303 : ! Initial guess for psi1
304 522 : hpsi0_ptr => p_psi0(:, idir)
305 : !
306 : !
307 522 : linres_control%converged = .FALSE.
308 522 : CALL linres_solver(p_env, qs_env, psi1_p(:, idir), hpsi0_ptr, psi0_order, output_unit, should_stop)
309 : !
310 : !
311 : ! print response functions
312 522 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
313 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
314 0 : ncubes = SIZE(list_cubes, 1)
315 0 : print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
316 0 : append_cube = section_get_lval(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%APPEND")
317 0 : my_pos = "REWIND"
318 0 : IF (append_cube) THEN
319 0 : my_pos = "APPEND"
320 : END IF
321 : !
322 :
323 0 : DO ispin = 1, nspins
324 : CALL qs_print_cubes(qs_env, psi1_p(ispin, idir), ncubes, list_cubes, &
325 : centers_set(ispin)%array, print_key, 'psi1_p', &
326 0 : idir=idir, ispin=ispin, file_position=my_pos)
327 : END DO ! ispin
328 : END IF ! print response functions
329 : !
330 : ! write restart file
331 696 : CALL linres_write_restart(qs_env, lr_section, psi1_p(:, idir), idir, "nmr_p")
332 : END DO ! idir
333 : !
334 : ! operator rxp
335 696 : DO idir = 1, 3
336 522 : IF (should_stop) EXIT
337 522 : IF (output_unit > 0) THEN
338 261 : WRITE (output_unit, "(T10,A)") "Response to the perturbation operator L_"//ACHAR(idir + 119)
339 : END IF
340 : !
341 : ! Initial guess for psi1
342 522 : hpsi0_ptr => rxp_psi0(:, idir)
343 : !
344 : !
345 522 : linres_control%converged = .FALSE.
346 522 : CALL linres_solver(p_env, qs_env, psi1_rxp(:, idir), hpsi0_ptr, psi0_order, output_unit, should_stop)
347 : !
348 : ! print response functions
349 522 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
350 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
351 0 : ncubes = SIZE(list_cubes, 1)
352 0 : print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
353 :
354 0 : DO ispin = 1, nspins
355 : CALL qs_print_cubes(qs_env, psi1_rxp(ispin, idir), ncubes, list_cubes, &
356 : centers_set(ispin)%array, print_key, 'psi1_rxp', &
357 0 : idir=idir, ispin=ispin, file_position=my_pos)
358 : END DO ! ispin
359 : END IF ! print response functions
360 : !
361 : ! write restart file
362 696 : CALL linres_write_restart(qs_env, lr_section, psi1_rxp(:, idir), idir, "nmr_rxp")
363 : END DO ! idir
364 174 : IF (.NOT. should_stop) current_env%all_pert_op_done = .TRUE.
365 : !
366 : ! operator D
367 174 : IF (current_env%full) THEN
368 142 : current_env%all_pert_op_done = .FALSE.
369 : !
370 : !
371 348 : DO ispin = 1, nspins
372 348 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
373 : END DO
374 : !
375 : ! The correction is state depedent a loop over the states is necessary
376 348 : max_nbr_center = MAXVAL(nbr_center(1:nspins))
377 348 : max_states = MAXVAL(nstates(1:nspins))
378 : !
379 0 : ALLOCATE (vecbuf_dklxp0(1, nao), fm_work_ii(nspins), fm_work_iii(nspins), &
380 1690 : dkl_vec_ii(max_states), dkl_vec_iii(max_states))
381 142 : vecbuf_dklxp0(1, nao) = 0.0_dp
382 : !
383 348 : DO ispin = 1, nspins
384 206 : nmo = nstates(ispin)
385 206 : mo_coeff => psi0_order(ispin)
386 206 : NULLIFY (tmp_fm_struct)
387 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
388 : ncol_global=nmo, para_env=para_env, &
389 206 : context=mo_coeff%matrix_struct%context)
390 :
391 206 : CALL cp_fm_create(fm_work_ii(ispin), tmp_fm_struct)
392 206 : CALL cp_fm_set_all(fm_work_ii(ispin), 0.0_dp)
393 206 : CALL cp_fm_create(fm_work_iii(ispin), tmp_fm_struct)
394 206 : CALL cp_fm_set_all(fm_work_iii(ispin), 0.0_dp)
395 348 : CALL cp_fm_struct_release(tmp_fm_struct)
396 : END DO ! ispin
397 : !
398 568 : DO idir = 1, 3
399 426 : IF (should_stop) EXIT
400 1044 : DO ispin = 1, nspins
401 1044 : CALL cp_fm_set_all(psi1_D(ispin, idir), 0.0_dp)
402 : END DO
403 426 : first_center = 0
404 426 : IF (linres_control%linres_restart) THEN
405 18 : CALL linres_read_restart(qs_env, lr_section, psi1_D(:, idir), idir, "nmr_dxp-", ind=first_center)
406 : END IF
407 426 : IF (first_center > 0) THEN
408 6 : IF (output_unit > 0) THEN
409 : WRITE (output_unit, "(T10,A,I4,A)")&
410 3 : & "Response to the perturbation operators (dk-dl)xp up to state ", &
411 6 : first_center, " have been read from restart"
412 : END IF
413 : END IF
414 : !
415 : ! here we run over the max number of states, then we need
416 : ! to be careful with overflow while doing uks calculations.
417 2886 : DO icenter = 1, max_nbr_center
418 : !
419 2460 : IF (should_stop) EXIT
420 2460 : IF (icenter > first_center) THEN
421 2436 : IF (output_unit > 0) THEN
422 : WRITE (output_unit, "(T10,A,I4,A)")&
423 1218 : & "Response to the perturbation operator (dk-dl)xp for -state- ", &
424 2436 : icenter, " in dir. "//ACHAR(idir + 119)
425 : END IF
426 : !
427 6306 : DO ispin = 1, nspins
428 3870 : nmo = nstates(ispin)
429 3870 : mo_coeff => psi0_order(ispin)
430 : !
431 : ! take care that no overflow can occur for uks
432 3870 : IF (icenter .GT. nbr_center(ispin)) THEN
433 : !
434 : ! set h1_psi0 and psi1 to zero to avoid problems in linres_scf
435 276 : CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
436 276 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
437 276 : CYCLE
438 : END IF
439 : !
440 30702 : dkl_vec_ii(:) = 0.0_dp
441 30702 : dkl_vec_iii(:) = 0.0_dp
442 : !
443 3594 : ist_true = statetrueindex(idir, icenter, ispin)
444 : !
445 : ! the initial guess is the previous set of psi1, just optimized
446 3594 : CALL set_vecp(idir, ii, iii)
447 14376 : dk(1:3) = centers_set(ispin)%array(1:3, ist_true)
448 : !
449 27612 : DO jcenter = 1, nbr_center(ispin)
450 96072 : dl(1:3) = centers_set(ispin)%array(1:3, jcenter)
451 24018 : dkl = pbc(dl, dk, cell)
452 52962 : DO j = center_list(ispin)%array(1, jcenter), center_list(ispin)%array(1, jcenter + 1) - 1
453 25350 : jstate = center_list(ispin)%array(2, j)
454 25350 : dkl_vec_ii(jstate) = dkl(ii)
455 49368 : dkl_vec_iii(jstate) = dkl(iii)
456 : END DO
457 : END DO
458 : !
459 : ! First term
460 : ! Rescale the ground state orbitals by (dk-dl)_ii
461 3594 : CALL cp_fm_to_fm(mo_coeff, fm_work_ii(ispin))
462 3594 : CALL cp_fm_column_scale(fm_work_ii(ispin), dkl_vec_ii(1:nmo))
463 : !
464 : ! Apply the p_iii operator
465 : ! fm_work_iii = -p_iii * (dk-dl)_ii * C0
466 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(iii)%matrix, fm_work_ii(ispin), &
467 3594 : fm_work_iii(ispin), ncol=nmo, alpha=-1.0_dp)
468 : !
469 : ! Copy in h1_psi0
470 : ! h1_psi0_i = fm_work_iii
471 3594 : CALL cp_fm_to_fm(fm_work_iii(ispin), h1_psi0(ispin))
472 : !
473 : ! Second term
474 : ! Rescale the ground state orbitals by (dk-dl)_iii
475 3594 : CALL cp_fm_to_fm(mo_coeff, fm_work_iii(ispin))
476 3594 : CALL cp_fm_column_scale(fm_work_iii(ispin), dkl_vec_iii(1:nmo))
477 : !
478 : ! Apply the p_ii operator
479 : ! fm_work_ii = -p_ii * (dk-dl)_iii * C0
480 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(ii)%matrix, fm_work_iii(ispin), &
481 3594 : fm_work_ii(ispin), ncol=nmo, alpha=-1.0_dp)
482 : !
483 : ! Copy in h1_psi0
484 : ! h1_psi0_i = fm_work_iii - fm_work_ii
485 : CALL cp_fm_scale_and_add(1.0_dp, h1_psi0(ispin),&
486 9624 : & -1.0_dp, fm_work_ii(ispin))
487 : END DO
488 :
489 : !
490 : ! Optimize the response wavefunctions
491 2436 : CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
492 : !
493 2436 : IF (output_unit > 0) THEN
494 : WRITE (output_unit, "(T10,A,/)")&
495 1218 : & "Store the psi1 vector for the calculation of the response current density "
496 : END IF
497 : !
498 6306 : DO ispin = 1, nspins
499 : !
500 : ! take care that no overflow can occur for uks
501 3870 : IF (icenter .GT. nbr_center(ispin)) CYCLE
502 : !
503 : ! need to reset those guys
504 3594 : ist_true = statetrueindex(idir, icenter, ispin)
505 7566 : DO i = center_list(ispin)%array(1, ist_true), center_list(ispin)%array(1, ist_true + 1) - 1
506 3972 : istate = center_list(ispin)%array(2, i)
507 : !
508 : ! the optimized wfns are copied in the fm
509 7566 : CALL cp_fm_to_fm(psi1(ispin), psi1_D(ispin, idir), 1, istate, istate)
510 : END DO
511 6306 : current_env%full_done(idir*ispin, icenter) = .TRUE.
512 : END DO ! ispin
513 : !
514 : ELSE
515 48 : DO ispin = 1, nspins
516 48 : current_env%full_done(idir*ispin, icenter) = .TRUE.
517 : END DO ! ispin
518 : END IF
519 2886 : CALL linres_write_restart(qs_env, lr_section, psi1_D(:, idir), idir, "nmr_dxp-", ind=icenter)
520 :
521 : END DO ! center
522 : !
523 : ! print response functions
524 426 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
525 142 : & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
526 0 : ncubes = SIZE(list_cubes, 1)
527 0 : print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
528 0 : DO ispin = 1, nspins
529 : CALL qs_print_cubes(qs_env, psi1_D(ispin, idir), &
530 : ncubes, list_cubes, centers_set(ispin)%array, print_key, 'psi1_D', &
531 0 : idir=idir, ispin=ispin, file_position=my_pos)
532 : END DO
533 : END IF ! print response functions
534 : !
535 : END DO ! idir
536 142 : IF (.NOT. should_stop) current_env%all_pert_op_done = .TRUE.
537 : !
538 : ! clean up
539 142 : CALL cp_fm_release(fm_work_ii)
540 142 : CALL cp_fm_release(fm_work_iii)
541 142 : DEALLOCATE (dkl_vec_ii, dkl_vec_iii, vecbuf_dklxp0)
542 :
543 : END IF
544 : !
545 : ! clean up
546 174 : IF (current_env%full) THEN
547 142 : CALL dbcsr_deallocate_matrix_set(op_p_ao)
548 142 : CALL cp_fm_release(psi1)
549 142 : CALL cp_fm_release(h1_psi0)
550 : END IF
551 : !
552 : CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
553 174 : & "PRINT%PROGRAM_RUN_INFO")
554 : !
555 174 : CALL timestop(handle)
556 : !
557 348 : END SUBROUTINE current_response
558 :
559 : ! **************************************************************************************************
560 :
561 : ! **************************************************************************************************
562 : !> \brief ...
563 : !> \param current_env ...
564 : !> \param qs_env ...
565 : ! **************************************************************************************************
566 174 : SUBROUTINE current_env_init(current_env, qs_env)
567 : !
568 : TYPE(current_env_type) :: current_env
569 : TYPE(qs_environment_type), POINTER :: qs_env
570 :
571 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_env_init'
572 :
573 : INTEGER :: handle, homo, i, iao, iatom, ibox, icenter, icount, idir, ii, ini, ir, is, ispin, &
574 : istate, istate2, istate_next, ix, iy, iz, j, jstate, k, max_nbr_center, max_states, n, &
575 : n_rep, nao, natom, nbr_box, ncubes, nmo, nspins, nstate, nstate_list(2), output_unit
576 174 : INTEGER, ALLOCATABLE, DIMENSION(:) :: buff, first_sgf, last_sgf
577 174 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
578 174 : INTEGER, DIMENSION(:), POINTER :: bounds, list, nbox, &
579 174 : selected_states_on_atom_list
580 : LOGICAL :: force_no_full, gapw, is0, &
581 : uniform_occupation
582 174 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: state_done
583 : REAL(dp) :: center(3), center2(3), dist, mdist, &
584 : r(3), rab(3)
585 174 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rbuff
586 174 : REAL(dp), DIMENSION(:), POINTER :: common_center
587 174 : REAL(dp), DIMENSION(:, :), POINTER :: center_array
588 174 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
589 : TYPE(cell_type), POINTER :: cell
590 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
591 : TYPE(cp_fm_type), POINTER :: mo_coeff
592 : TYPE(cp_logger_type), POINTER :: logger
593 : TYPE(dft_control_type), POINTER :: dft_control
594 174 : TYPE(jrho_atom_type), DIMENSION(:), POINTER :: jrho1_atom_set
595 : TYPE(linres_control_type), POINTER :: linres_control
596 : TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control
597 174 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
598 : TYPE(mp_para_env_type), POINTER :: para_env
599 174 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
600 174 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
601 : TYPE(pw_env_type), POINTER :: pw_env
602 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
603 174 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
604 174 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
605 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
606 : TYPE(qs_matrix_pools_type), POINTER :: mpools
607 : TYPE(scf_control_type), POINTER :: scf_control
608 : TYPE(section_vals_type), POINTER :: current_section, lr_section
609 :
610 174 : CALL timeset(routineN, handle)
611 :
612 174 : NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, linres_control, scf_control, &
613 174 : logger, mos, mpools, current_section, particle_set, mo_coeff, &
614 174 : auxbas_pw_pool, pw_env, jrho1_atom_set, common_center, tmp_fm_struct, &
615 174 : para_env, qs_loc_env, localized_wfn_control, rho_g, rho_r)
616 : !
617 : !
618 : CALL get_qs_env(qs_env=qs_env, &
619 : atomic_kind_set=atomic_kind_set, &
620 : qs_kind_set=qs_kind_set, &
621 : cell=cell, &
622 : dft_control=dft_control, &
623 : linres_control=linres_control, &
624 : mos=mos, &
625 : mpools=mpools, &
626 : particle_set=particle_set, &
627 : pw_env=pw_env, &
628 : scf_control=scf_control, &
629 174 : para_env=para_env)
630 : !
631 174 : gapw = dft_control%qs_control%gapw
632 174 : nspins = dft_control%nspins
633 174 : natom = SIZE(particle_set, 1)
634 174 : CALL get_mo_set(mo_set=mos(1), nao=nao)
635 : !
636 174 : max_states = 0
637 424 : DO ispin = 1, nspins
638 250 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
639 424 : max_states = MAX(max_states, nmo)
640 : END DO
641 : !
642 : !
643 : !
644 : ! some checks
645 424 : DO ispin = 1, nspins
646 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, homo=homo, &
647 250 : uniform_occupation=uniform_occupation)
648 : !
649 : ! check that homo=nmo
650 250 : IF (nmo .NE. homo) CPABORT("nmo != homo")
651 : !
652 : ! check that the nbr of localized states is equal to the nbr of states
653 : ! nmoloc = SIZE(linres_control%localized_wfn_control%centers_set(ispin)%array,2)
654 : ! IF (nmoloc.NE.nmo) CALL cp_abort(__LOCATION__,&
655 : ! "The nbr of localized functions "//&
656 : ! & "is not equal to the nbr of states.")
657 : !
658 : ! check that the occupation is uniform
659 674 : IF (.NOT. uniform_occupation) CPABORT("nmo != homo")
660 : END DO
661 : !
662 : !
663 174 : logger => cp_get_default_logger()
664 174 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
665 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
666 174 : extension=".linresLog")
667 :
668 174 : IF (output_unit > 0) THEN
669 87 : WRITE (output_unit, "(/,T20,A,/)") "*** Start current Calculation ***"
670 87 : WRITE (output_unit, "(T10,A,/)") "Inizialization of the current environment"
671 : END IF
672 :
673 174 : CALL current_env_cleanup(current_env)
674 174 : current_env%gauge_init = .FALSE.
675 4698 : current_env%chi_tensor(:, :, :) = 0.0_dp
676 4698 : current_env%chi_tensor_loc(:, :, :) = 0.0_dp
677 174 : current_env%nao = nao
678 174 : current_env%full = .TRUE.
679 174 : current_env%do_selected_states = .FALSE.
680 : !
681 : ! If current_density or full_nmr different allocations are required
682 174 : current_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%CURRENT")
683 : !
684 : ! Select the gauge
685 174 : CALL section_vals_val_get(current_section, "GAUGE", i_val=current_env%gauge)
686 178 : SELECT CASE (current_env%gauge)
687 : CASE (current_gauge_r)
688 4 : current_env%gauge_name = "R"
689 : CASE (current_gauge_r_and_step_func)
690 142 : current_env%gauge_name = "R_AND_STEP_FUNCTION"
691 : CASE (current_gauge_atom)
692 28 : current_env%gauge_name = "ATOM"
693 : CASE DEFAULT
694 174 : CPABORT("Unknown gauge, try again...")
695 : END SELECT
696 : !
697 : ! maximal radius to build the atom gauge
698 : CALL section_vals_val_get(current_section, "GAUGE_ATOM_RADIUS", &
699 174 : r_val=current_env%gauge_atom_radius)
700 : ! use old gauge atom
701 : CALL section_vals_val_get(current_section, "USE_OLD_GAUGE_ATOM", &
702 174 : l_val=current_env%use_old_gauge_atom)
703 : ! chi for pbc
704 174 : CALL section_vals_val_get(current_section, "CHI_PBC", l_val=current_env%chi_pbc)
705 : !
706 : ! use old gauge atom
707 : CALL section_vals_val_get(current_section, "USE_OLD_GAUGE_ATOM", &
708 174 : l_val=current_env%use_old_gauge_atom)
709 : !
710 : ! chi for pbc
711 174 : CALL section_vals_val_get(current_section, "CHI_PBC", l_val=current_env%chi_pbc)
712 : !
713 : ! which center for the orbitals shall we use
714 174 : CALL section_vals_val_get(current_section, "ORBITAL_CENTER", i_val=current_env%orb_center)
715 300 : SELECT CASE (current_env%orb_center)
716 : CASE (current_orb_center_wannier)
717 : !
718 126 : current_env%orb_center_name = "WANNIER"
719 : CASE (current_orb_center_common)
720 : !
721 32 : current_env%orb_center_name = "COMMON"
722 32 : current_env%full = .FALSE.
723 : !
724 : ! Is there a user specified common_center?
725 32 : CALL section_vals_val_get(current_section, "COMMON_CENTER", r_vals=common_center)
726 : CASE (current_orb_center_atom)
727 : !
728 10 : current_env%orb_center_name = "ATOM"
729 : CASE (current_orb_center_box)
730 : !
731 6 : current_env%orb_center_name = "BOX"
732 : !
733 : ! Is there a user specified nbox?
734 6 : CALL section_vals_val_get(current_section, "NBOX", i_vals=nbox)
735 : CASE DEFAULT
736 174 : CPABORT("Unknown orbital center, try again...")
737 : END SELECT
738 :
739 174 : CALL section_vals_val_get(current_section, "FORCE_NO_FULL", l_val=force_no_full)
740 174 : IF (force_no_full) current_env%full = .FALSE.
741 : !
742 : ! Check if restat also psi0 should be restarted
743 : !IF(current_env%restart_current .AND. scf_control%density_guess/=restart_guess)THEN
744 : ! CPABORT("restart_nmr requires density_guess=restart")
745 : !ENDIF
746 : !
747 : ! check that the psi0 are localized and you have all the centers
748 174 : CPASSERT(linres_control%localized_psi0)
749 174 : IF (output_unit > 0) THEN
750 : WRITE (output_unit, '(A)') &
751 87 : ' To get CURRENT parameters within PBC you need localized zero order orbitals '
752 : END IF
753 174 : qs_loc_env => linres_control%qs_loc_env
754 174 : CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
755 : !
756 : !
757 : ALLOCATE (current_env%centers_set(nspins), current_env%center_list(nspins), &
758 1718 : state_list(max_states, nspins))
759 2196 : state_list(:, :) = HUGE(0)
760 522 : nstate_list(:) = HUGE(0)
761 : !
762 : !
763 : !
764 : ! if qmmm is selected change the definition of the center of the box 0 -> L/2
765 174 : IF (current_env%do_qmmm) THEN
766 12 : DO ispin = 1, nspins
767 42 : DO istate = 1, SIZE(localized_wfn_control%centers_set(ispin)%array, 2)
768 : ! just to be sure...
769 30 : r(:) = pbc(localized_wfn_control%centers_set(ispin)%array(:, istate), cell)
770 30 : IF (r(1) .LT. 0.0_dp) THEN
771 : localized_wfn_control%centers_set(ispin)%array(1, istate) = &
772 6 : r(1) + cell%hmat(1, 1)
773 : END IF
774 30 : IF (r(2) .LT. 0.0_dp) THEN
775 : localized_wfn_control%centers_set(ispin)%array(2, istate) = &
776 24 : r(2) + cell%hmat(2, 2)
777 : END IF
778 36 : IF (r(3) .LT. 0.0_dp) THEN
779 : localized_wfn_control%centers_set(ispin)%array(3, istate) = &
780 6 : r(3) + cell%hmat(3, 3)
781 : END IF
782 : END DO
783 : END DO
784 : END IF
785 : !
786 : !
787 : !
788 : ! if the user has requested to compute the response for a subset of the states
789 : ! we collect them here. it requies the states to be localized!
790 : CALL section_vals_val_get(current_section, "SELECTED_STATES_ATOM_RADIUS", &
791 174 : r_val=current_env%selected_states_atom_radius)
792 174 : CALL section_vals_val_get(current_section, "SELECTED_STATES_ON_ATOM_LIST", n_rep_val=n_rep)
793 : !
794 174 : current_env%do_selected_states = n_rep .GT. 0
795 : !
796 : ! for the moment selected states doesnt work with the preconditioner FULL_ALL
797 174 : IF (linres_control%preconditioner_type .EQ. ot_precond_full_all .AND. current_env%do_selected_states) &
798 0 : CPABORT("Selected states doesnt work with the preconditioner FULL_ALL")
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)) .LE. current_env%selected_states_atom_radius
827 1456 : IF ((DOT_PRODUCT(r, r)) .LE. (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 158 : 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 .LT. 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) .EQ. 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 96 : 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 .LT. 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) .EQ. 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 .EQ. current_orb_center_wannier .OR. &
1051 : current_env%orb_center .EQ. current_orb_center_atom .OR. &
1052 : current_env%orb_center .EQ. 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) .LE. 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) .LE. 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) .LE. 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 5618 : 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 .LT. mdist) THEN
1124 6790 : mdist = dist
1125 6790 : 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 .EQ. 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 .EQ. current_orb_center_common) THEN
1162 16 : WRITE (output_unit, "(T2,A,T50,3F10.6)") "CURRENT| Common center", common_center(1:3)
1163 71 : ELSEIF (current_env%orb_center .EQ. current_orb_center_box) THEN
1164 3 : 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 174 : 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
|
