(git:e7e05ae)
qs_linres_current_utils.F
Go to the documentation of this file.
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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
140  INTEGER, DIMENSION(:), POINTER :: list_cubes, row_blk_sizes
141  INTEGER, DIMENSION(:, :, :), POINTER :: statetrueindex
142  LOGICAL :: append_cube, should_stop
143  REAL(dp) :: dk(3), dkl(3), dl(3)
144  REAL(dp), ALLOCATABLE, DIMENSION(:) :: dkl_vec_ii, dkl_vec_iii
145  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vecbuf_dklxp0
146  TYPE(cell_type), POINTER :: cell
147  TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
148  TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
149  TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
150  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_work_ii, fm_work_iii, h1_psi0, psi1
151  TYPE(cp_fm_type), DIMENSION(:), POINTER :: hpsi0_ptr, psi0_order
152  TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: p_psi0, psi1_d, psi1_p, psi1_rxp, &
153  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  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  POINTER :: sab_orb
163  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
164  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  CALL timeset(routinen, handle)
169  !
170  NULLIFY (cell, dft_control, linres_control, lr_section, current_section, &
171  logger, mpools, mo_coeff, para_env, &
172  tmp_fm_struct, &
173  list_cubes, statetrueindex, centers_set, center_list, psi1_p, psi1_rxp, psi1_d, &
174  p_psi0, rxp_psi0, psi0_order, op_p_ao, sab_orb, particle_set)
175 
176  logger => cp_get_default_logger()
177  lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
178  current_section => section_vals_get_subs_vals(qs_env%input, &
179  "PROPERTIES%LINRES%CURRENT")
180 
181  output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
182  extension=".linresLog")
183  IF (output_unit > 0) THEN
184  WRITE (unit=output_unit, fmt="(T10,A,/)") &
185  "*** 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  para_env=para_env)
197 
198  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  psi0_order=psi0_order)
214  !
215  ! allocate the vectors
216  IF (current_env%full) THEN
217  ALLOCATE (psi1(nspins), h1_psi0(nspins))
218  DO ispin = 1, nspins
219  mo_coeff => psi0_order(ispin)
220  NULLIFY (tmp_fm_struct)
221  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
222  ncol_global=nstates(ispin), &
223  context=mo_coeff%matrix_struct%context)
224  CALL cp_fm_create(psi1(ispin), tmp_fm_struct)
225  CALL cp_fm_create(h1_psi0(ispin), tmp_fm_struct)
226  CALL cp_fm_struct_release(tmp_fm_struct)
227  END DO
228  !
229  ! prepare for allocation
230  natom = SIZE(particle_set, 1)
231  ALLOCATE (first_sgf(natom))
232  ALLOCATE (last_sgf(natom))
233  CALL get_particle_set(particle_set, qs_kind_set, &
234  first_sgf=first_sgf, &
235  last_sgf=last_sgf)
236  ALLOCATE (row_blk_sizes(natom))
237  CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
238  DEALLOCATE (first_sgf, last_sgf)
239  !
240  ! rebuild the linear momentum matrices
241  CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
242  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  nze=0, mutable_work=.true.)
248  CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
249  !
250  !
251  DEALLOCATE (row_blk_sizes)
252  !
253  !
254  CALL dbcsr_set(op_p_ao(1)%matrix, 0.0_dp)
255  DO idir = 2, 3
256  ALLOCATE (op_p_ao(idir)%matrix)
257  CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
258  "current_env%op_p_ao"//"-"//trim(adjustl(cp_to_string(idir))))
259  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  CALL build_lin_mom_matrix(qs_env, op_p_ao)
264  !
265  END IF
266  !
267  ! set response to zero
268  !
269  DO idir = 1, 3
270  DO ispin = 1, nspins
271  CALL cp_fm_set_all(psi1_p(ispin, idir), 0.0_dp)
272  CALL cp_fm_set_all(psi1_rxp(ispin, idir), 0.0_dp)
273  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  first_center = 0
282  IF (linres_control%linres_restart) THEN
283  DO idir = 1, 3
284  ! operator p
285  CALL linres_read_restart(qs_env, lr_section, psi1_p(:, idir), idir, "nmr_p")
286  ! operator rxp
287  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  should_stop = .false.
295  current_env%all_pert_op_done = .false.
296  ! operator p
297  DO idir = 1, 3
298  IF (should_stop) EXIT
299  IF (output_unit > 0) THEN
300  WRITE (output_unit, "(T10,A)") "Response to the perturbation operator P_"//achar(idir + 119)
301  END IF
302  !
303  ! Initial guess for psi1
304  hpsi0_ptr => p_psi0(:, idir)
305  !
306  !
307  linres_control%converged = .false.
308  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  IF (btest(cp_print_key_should_output(logger%iter_info, current_section,&
313  & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
314  ncubes = SIZE(list_cubes, 1)
315  print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
316  append_cube = section_get_lval(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%APPEND")
317  my_pos = "REWIND"
318  IF (append_cube) THEN
319  my_pos = "APPEND"
320  END IF
321  !
322 
323  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  idir=idir, ispin=ispin, file_position=my_pos)
327  END DO ! ispin
328  END IF ! print response functions
329  !
330  ! write restart file
331  CALL linres_write_restart(qs_env, lr_section, psi1_p(:, idir), idir, "nmr_p")
332  END DO ! idir
333  !
334  ! operator rxp
335  DO idir = 1, 3
336  IF (should_stop) EXIT
337  IF (output_unit > 0) THEN
338  WRITE (output_unit, "(T10,A)") "Response to the perturbation operator L_"//achar(idir + 119)
339  END IF
340  !
341  ! Initial guess for psi1
342  hpsi0_ptr => rxp_psi0(:, idir)
343  !
344  !
345  linres_control%converged = .false.
346  CALL linres_solver(p_env, qs_env, psi1_rxp(:, idir), hpsi0_ptr, psi0_order, output_unit, should_stop)
347  !
348  ! print response functions
349  IF (btest(cp_print_key_should_output(logger%iter_info, current_section,&
350  & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
351  ncubes = SIZE(list_cubes, 1)
352  print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
353 
354  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  idir=idir, ispin=ispin, file_position=my_pos)
358  END DO ! ispin
359  END IF ! print response functions
360  !
361  ! write restart file
362  CALL linres_write_restart(qs_env, lr_section, psi1_rxp(:, idir), idir, "nmr_rxp")
363  END DO ! idir
364  IF (.NOT. should_stop) current_env%all_pert_op_done = .true.
365  !
366  ! operator D
367  IF (current_env%full) THEN
368  current_env%all_pert_op_done = .false.
369  !
370  !
371  DO ispin = 1, nspins
372  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  max_nbr_center = maxval(nbr_center(1:nspins))
377  max_states = maxval(nstates(1:nspins))
378  !
379  ALLOCATE (vecbuf_dklxp0(1, nao), fm_work_ii(nspins), fm_work_iii(nspins), &
380  dkl_vec_ii(max_states), dkl_vec_iii(max_states))
381  vecbuf_dklxp0(1, nao) = 0.0_dp
382  !
383  DO ispin = 1, nspins
384  nmo = nstates(ispin)
385  mo_coeff => psi0_order(ispin)
386  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  context=mo_coeff%matrix_struct%context)
390 
391  CALL cp_fm_create(fm_work_ii(ispin), tmp_fm_struct)
392  CALL cp_fm_set_all(fm_work_ii(ispin), 0.0_dp)
393  CALL cp_fm_create(fm_work_iii(ispin), tmp_fm_struct)
394  CALL cp_fm_set_all(fm_work_iii(ispin), 0.0_dp)
395  CALL cp_fm_struct_release(tmp_fm_struct)
396  END DO ! ispin
397  !
398  DO idir = 1, 3
399  IF (should_stop) EXIT
400  DO ispin = 1, nspins
401  CALL cp_fm_set_all(psi1_d(ispin, idir), 0.0_dp)
402  END DO
403  first_center = 0
404  IF (linres_control%linres_restart) THEN
405  CALL linres_read_restart(qs_env, lr_section, psi1_d(:, idir), idir, "nmr_dxp-", ind=first_center)
406  END IF
407  IF (first_center > 0) THEN
408  IF (output_unit > 0) THEN
409  WRITE (output_unit, "(T10,A,I4,A)")&
410  & "Response to the perturbation operators (dk-dl)xp up to state ", &
411  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  DO icenter = 1, max_nbr_center
418  !
419  IF (should_stop) EXIT
420  IF (icenter > first_center) THEN
421  IF (output_unit > 0) THEN
422  WRITE (output_unit, "(T10,A,I4,A)")&
423  & operator"Response to the perturbation (dk-dl)xp for -state- ", &
424  icenter, " in dir. "//achar(idir + 119)
425  END IF
426  !
427  DO ispin = 1, nspins
428  nmo = nstates(ispin)
429  mo_coeff => psi0_order(ispin)
430  !
431  ! take care that no overflow can occur for uks
432  IF (icenter .GT. nbr_center(ispin)) THEN
433  !
434  ! set h1_psi0 and psi1 to zero to avoid problems in linres_scf
435  CALL cp_fm_set_all(h1_psi0(ispin), 0.0_dp)
436  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
437  cycle
438  END IF
439  !
440  dkl_vec_ii(:) = 0.0_dp
441  dkl_vec_iii(:) = 0.0_dp
442  !
443  ist_true = statetrueindex(idir, icenter, ispin)
444  !
445  ! the initial guess is the previous set of psi1, just optimized
446  CALL set_vecp(idir, ii, iii)
447  dk(1:3) = centers_set(ispin)%array(1:3, ist_true)
448  !
449  DO jcenter = 1, nbr_center(ispin)
450  dl(1:3) = centers_set(ispin)%array(1:3, jcenter)
451  dkl = pbc(dl, dk, cell)
452  DO j = center_list(ispin)%array(1, jcenter), center_list(ispin)%array(1, jcenter + 1) - 1
453  jstate = center_list(ispin)%array(2, j)
454  dkl_vec_ii(jstate) = dkl(ii)
455  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  CALL cp_fm_to_fm(mo_coeff, fm_work_ii(ispin))
462  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  fm_work_iii(ispin), ncol=nmo, alpha=-1.0_dp)
468  !
469  ! Copy in h1_psi0
470  ! h1_psi0_i = fm_work_iii
471  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  CALL cp_fm_to_fm(mo_coeff, fm_work_iii(ispin))
476  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  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  & -1.0_dp, fm_work_ii(ispin))
487  END DO
488 
489  !
490  ! Optimize the response wavefunctions
491  CALL linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, output_unit, should_stop)
492  !
493  IF (output_unit > 0) THEN
494  WRITE (output_unit, "(T10,A,/)")&
495  & "Store the psi1 vector for the calculation of the response current density "
496  END IF
497  !
498  DO ispin = 1, nspins
499  !
500  ! take care that no overflow can occur for uks
501  IF (icenter .GT. nbr_center(ispin)) cycle
502  !
503  ! need to reset those guys
504  ist_true = statetrueindex(idir, icenter, ispin)
505  DO i = center_list(ispin)%array(1, ist_true), center_list(ispin)%array(1, ist_true + 1) - 1
506  istate = center_list(ispin)%array(2, i)
507  !
508  ! the optimized wfns are copied in the fm
509  CALL cp_fm_to_fm(psi1(ispin), psi1_d(ispin, idir), 1, istate, istate)
510  END DO
511  current_env%full_done(idir*ispin, icenter) = .true.
512  END DO ! ispin
513  !
514  ELSE
515  DO ispin = 1, nspins
516  current_env%full_done(idir*ispin, icenter) = .true.
517  END DO ! ispin
518  END IF
519  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  IF (btest(cp_print_key_should_output(logger%iter_info, current_section,&
525  & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
526  ncubes = SIZE(list_cubes, 1)
527  print_key => section_vals_get_subs_vals(current_section, "PRINT%RESPONSE_FUNCTION_CUBES")
528  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  idir=idir, ispin=ispin, file_position=my_pos)
532  END DO
533  END IF ! print response functions
534  !
535  END DO ! idir
536  IF (.NOT. should_stop) current_env%all_pert_op_done = .true.
537  !
538  ! clean up
539  CALL cp_fm_release(fm_work_ii)
540  CALL cp_fm_release(fm_work_iii)
541  DEALLOCATE (dkl_vec_ii, dkl_vec_iii, vecbuf_dklxp0)
542 
543  END IF
544  !
545  ! clean up
546  IF (current_env%full) THEN
547  CALL dbcsr_deallocate_matrix_set(op_p_ao)
548  CALL cp_fm_release(psi1)
549  CALL cp_fm_release(h1_psi0)
550  END IF
551  !
552  CALL cp_print_key_finished_output(output_unit, logger, lr_section,&
553  & "PRINT%PROGRAM_RUN_INFO")
554  !
555  CALL timestop(handle)
556  !
557  END SUBROUTINE current_response
558 
559 ! **************************************************************************************************
560 
561 ! **************************************************************************************************
562 !> \brief ...
563 !> \param current_env ...
564 !> \param qs_env ...
565 ! **************************************************************************************************
566  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: buff, first_sgf, last_sgf
577  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: state_list
578  INTEGER, DIMENSION(:), POINTER :: bounds, list, nbox, &
579  selected_states_on_atom_list
580  LOGICAL :: force_no_full, gapw, is0, &
581  uniform_occupation
582  LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: state_done
583  REAL(dp) :: center(3), center2(3), dist, mdist, &
584  r(3), rab(3)
585  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rbuff
586  REAL(dp), DIMENSION(:), POINTER :: common_center
587  REAL(dp), DIMENSION(:, :), POINTER :: center_array
588  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  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  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
598  TYPE(mp_para_env_type), POINTER :: para_env
599  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
600  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  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
604  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  CALL timeset(routinen, handle)
611 
612  NULLIFY (atomic_kind_set, qs_kind_set, cell, dft_control, linres_control, scf_control, &
613  logger, mos, mpools, current_section, particle_set, mo_coeff, &
614  auxbas_pw_pool, pw_env, jrho1_atom_set, common_center, tmp_fm_struct, &
615  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  para_env=para_env)
630  !
631  gapw = dft_control%qs_control%gapw
632  nspins = dft_control%nspins
633  natom = SIZE(particle_set, 1)
634  CALL get_mo_set(mo_set=mos(1), nao=nao)
635  !
636  max_states = 0
637  DO ispin = 1, nspins
638  CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
639  max_states = max(max_states, nmo)
640  END DO
641  !
642  !
643  !
644  ! some checks
645  DO ispin = 1, nspins
646  CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, homo=homo, &
647  uniform_occupation=uniform_occupation)
648  !
649  ! check that homo=nmo
650  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  IF (.NOT. uniform_occupation) cpabort("nmo != homo")
660  END DO
661  !
662  !
663  logger => cp_get_default_logger()
664  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  extension=".linresLog")
667 
668  IF (output_unit > 0) THEN
669  WRITE (output_unit, "(/,T20,A,/)") "*** Start current Calculation ***"
670  WRITE (output_unit, "(T10,A,/)") "Inizialization of the current environment"
671  END IF
672 
673  CALL current_env_cleanup(current_env)
674  current_env%gauge_init = .false.
675  current_env%chi_tensor(:, :, :) = 0.0_dp
676  current_env%chi_tensor_loc(:, :, :) = 0.0_dp
677  current_env%nao = nao
678  current_env%full = .true.
679  current_env%do_selected_states = .false.
680  !
681  ! If current_density or full_nmr different allocations are required
682  current_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%CURRENT")
683  !
684  ! Select the gauge
685  CALL section_vals_val_get(current_section, "GAUGE", i_val=current_env%gauge)
686  SELECT CASE (current_env%gauge)
687  CASE (current_gauge_r)
688  current_env%gauge_name = "R"
689  CASE (current_gauge_r_and_step_func)
690  current_env%gauge_name = "R_AND_STEP_FUNCTION"
691  CASE (current_gauge_atom)
692  current_env%gauge_name = "ATOM"
693  CASE DEFAULT
694  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  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  l_val=current_env%use_old_gauge_atom)
703  ! chi for pbc
704  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  l_val=current_env%use_old_gauge_atom)
709  !
710  ! chi for pbc
711  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  CALL section_vals_val_get(current_section, "ORBITAL_CENTER", i_val=current_env%orb_center)
715  SELECT CASE (current_env%orb_center)
716  CASE (current_orb_center_wannier)
717  !
718  current_env%orb_center_name = "WANNIER"
719  CASE (current_orb_center_common)
720  !
721  current_env%orb_center_name = "COMMON"
722  current_env%full = .false.
723  !
724  ! Is there a user specified common_center?
725  CALL section_vals_val_get(current_section, "COMMON_CENTER", r_vals=common_center)
726  CASE (current_orb_center_atom)
727  !
728  current_env%orb_center_name = "ATOM"
729  CASE (current_orb_center_box)
730  !
731  current_env%orb_center_name = "BOX"
732  !
733  ! Is there a user specified nbox?
734  CALL section_vals_val_get(current_section, "NBOX", i_vals=nbox)
735  CASE DEFAULT
736  cpabort("Unknown orbital center, try again...")
737  END SELECT
738 
739  CALL section_vals_val_get(current_section, "FORCE_NO_FULL", l_val=force_no_full)
740  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  cpassert(linres_control%localized_psi0)
749  IF (output_unit > 0) THEN
750  WRITE (output_unit, '(A)') &
751  ' To get CURRENT parameters within PBC you need localized zero order orbitals '
752  END IF
753  qs_loc_env => linres_control%qs_loc_env
754  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  state_list(max_states, nspins))
759  state_list(:, :) = huge(0)
760  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  IF (current_env%do_qmmm) THEN
766  DO ispin = 1, nspins
767  DO istate = 1, SIZE(localized_wfn_control%centers_set(ispin)%array, 2)
768  ! just to be sure...
769  r(:) = pbc(localized_wfn_control%centers_set(ispin)%array(:, istate), cell)
770  IF (r(1) .LT. 0.0_dp) THEN
771  localized_wfn_control%centers_set(ispin)%array(1, istate) = &
772  r(1) + cell%hmat(1, 1)
773  END IF
774  IF (r(2) .LT. 0.0_dp) THEN
775  localized_wfn_control%centers_set(ispin)%array(2, istate) = &
776  r(2) + cell%hmat(2, 2)
777  END IF
778  IF (r(3) .LT. 0.0_dp) THEN
779  localized_wfn_control%centers_set(ispin)%array(3, istate) = &
780  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  r_val=current_env%selected_states_atom_radius)
792  CALL section_vals_val_get(current_section, "SELECTED_STATES_ON_ATOM_LIST", n_rep_val=n_rep)
793  !
794  current_env%do_selected_states = n_rep .GT. 0
795  !
796  ! for the moment selected states doesnt work with the preconditioner FULL_ALL
797  IF (linres_control%preconditioner_type .EQ. ot_precond_full_all .AND. current_env%do_selected_states) &
798  cpabort("Selected states doesnt work with the preconditioner FULL_ALL")
799  !
800  !
801  NULLIFY (current_env%selected_states_on_atom_list)
802  n = 0
803  DO ir = 1, n_rep
804  NULLIFY (list)
805  CALL section_vals_val_get(current_section, "SELECTED_STATES_ON_ATOM_LIST", &
806  i_rep_val=ir, i_vals=list)
807  IF (ASSOCIATED(list)) THEN
808  CALL reallocate(current_env%selected_states_on_atom_list, 1, n + SIZE(list))
809  DO ini = 1, SIZE(list)
810  current_env%selected_states_on_atom_list(ini + n) = list(ini)
811  END DO
812  n = n + SIZE(list)
813  END IF
814  END DO
815  !
816  ! build the subset
817  IF (current_env%do_selected_states) THEN
818  selected_states_on_atom_list => current_env%selected_states_on_atom_list
819  DO ispin = 1, nspins
820  center_array => localized_wfn_control%centers_set(ispin)%array
821  nstate = 0
822  DO istate = 1, SIZE(center_array, 2)
823  DO i = 1, SIZE(selected_states_on_atom_list, 1)
824  iatom = selected_states_on_atom_list(i)
825  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  IF ((dot_product(r, r)) .LE. (current_env%selected_states_atom_radius &
828  *current_env%selected_states_atom_radius)) &
829  THEN
830  !
831  ! add the state to the list
832  nstate = nstate + 1
833  state_list(nstate, ispin) = istate
834  EXIT
835  END IF
836  END DO
837  END DO
838  nstate_list(ispin) = nstate
839  END DO
840  ELSE
841  DO ispin = 1, nspins
842  center_array => localized_wfn_control%centers_set(ispin)%array
843  nstate = 0
844  DO istate = 1, SIZE(center_array, 2)
845  nstate = nstate + 1
846  state_list(nstate, ispin) = istate
847  END DO
848  nstate_list(ispin) = nstate
849  END DO
850  END IF
851  !
852  !
853  !
854  ! clustering the states
855  DO ispin = 1, nspins
856  nstate = nstate_list(ispin)
857  current_env%nstates(ispin) = nstate
858  !
859  ALLOCATE (current_env%center_list(ispin)%array(2, nstate + 1), &
860  current_env%centers_set(ispin)%array(3, nstate))
861  current_env%center_list(ispin)%array(:, :) = huge(0)
862  current_env%centers_set(ispin)%array(:, :) = huge(0.0_dp)
863  !
864  center_array => localized_wfn_control%centers_set(ispin)%array
865  !
866  ! point to the psi0 centers
867  SELECT CASE (current_env%orb_center)
868  CASE (current_orb_center_wannier)
869  !
870  ! use the wannier center as -center-
871  current_env%nbr_center(ispin) = nstate
872  DO is = 1, nstate
873  istate = state_list(is, ispin)
874  current_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
875  current_env%center_list(ispin)%array(1, is) = is
876  current_env%center_list(ispin)%array(2, is) = istate
877  END DO
878  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  current_env%centers_set(ispin)%array(:, 1) = common_center(:)
884  current_env%nbr_center(ispin) = 1
885  current_env%center_list(ispin)%array(1, 1) = 1
886  current_env%center_list(ispin)%array(1, 2) = nstate + 1
887  DO is = 1, nstate
888  istate = state_list(is, ispin)
889  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  ALLOCATE (buff(nstate_list(ispin)))
896  buff(:) = 0
897  !
898  DO is = 1, nstate
899  istate = state_list(is, ispin)
900  mdist = huge(0.0_dp)
901  DO iatom = 1, natom
902  r = pbc(particle_set(iatom)%r(:), cell)
903  rab = pbc(r, center_array(1:3, istate), cell)
904  dist = sqrt(rab(1)**2 + rab(2)**2 + rab(3)**2)
905  IF (dist .LT. mdist) THEN
906  buff(is) = iatom
907  mdist = dist
908  END IF
909  END DO
910  END DO
911  !
912  i = 0
913  ii = 1
914  current_env%center_list(ispin)%array(1, 1) = 1
915  DO iatom = 1, natom
916  j = 0
917  is0 = .true.
918  DO is = 1, nstate
919  istate = state_list(is, ispin)
920  IF (buff(is) .EQ. iatom) THEN
921  j = j + 1
922  i = i + 1
923  is0 = .false.
924  current_env%center_list(ispin)%array(2, i) = istate
925  END IF
926  END DO
927  IF (.NOT. is0) THEN
928  IF (output_unit > 0) THEN
929  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  current_env%center_list(ispin)%array(1, ii) + j
933  current_env%centers_set(ispin)%array(:, ii) = &
934  pbc(particle_set(iatom)%r, cell)
935  ii = ii + 1
936  END IF
937  END DO
938  current_env%nbr_center(ispin) = ii - 1
939  !
940  DEALLOCATE (buff)
941  CASE (current_orb_center_box)
942  !
943  ! use boxes as -center-
944  nbr_box = nbox(1)*nbox(2)*nbox(3)
945  ALLOCATE (rbuff(3, nbr_box), buff(nstate))
946  rbuff(:, :) = huge(0.0_dp)
947  buff(:) = 0
948  !
949  ibox = 1
950  DO iz = 1, nbox(3)
951  DO iy = 1, nbox(2)
952  DO ix = 1, nbox(1)
953  rbuff(1, ibox) = cell%hmat(1, 1)*((real(ix, dp) - 0.5_dp)/real(nbox(1), dp) - 0.5_dp)
954  rbuff(2, ibox) = cell%hmat(2, 2)*((real(iy, dp) - 0.5_dp)/real(nbox(2), dp) - 0.5_dp)
955  rbuff(3, ibox) = cell%hmat(3, 3)*((real(iz, dp) - 0.5_dp)/real(nbox(3), dp) - 0.5_dp)
956  ibox = ibox + 1
957  END DO
958  END DO
959  END DO
960  !
961  DO is = 1, nstate
962  istate = state_list(is, ispin)
963  mdist = huge(0.0_dp)
964  DO ibox = 1, nbr_box
965  rab(:) = pbc(rbuff(:, ibox), center_array(1:3, istate), cell)
966  dist = sqrt(rab(1)**2 + rab(2)**2 + rab(3)**2)
967  IF (dist .LT. mdist) THEN
968  buff(is) = ibox
969  mdist = dist
970  END IF
971  END DO
972  END DO
973  !
974  i = 0
975  ii = 1
976  current_env%center_list(ispin)%array(1, 1) = 1
977  DO ibox = 1, nbr_box
978  j = 0
979  is0 = .true.
980  DO is = 1, nstate
981  istate = state_list(is, ispin)
982  IF (buff(is) .EQ. ibox) THEN
983  j = j + 1
984  i = i + 1
985  is0 = .false.
986  current_env%center_list(ispin)%array(2, i) = istate
987  END IF
988  END DO
989  IF (.NOT. is0) THEN
990  IF (output_unit > 0) THEN
991  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  current_env%center_list(ispin)%array(1, ii) + j
995  current_env%centers_set(ispin)%array(:, ii) = rbuff(:, ibox)
996  ii = ii + 1
997  END IF
998  END DO
999  current_env%nbr_center(ispin) = ii - 1
1000  !
1001  DEALLOCATE (buff, rbuff)
1002  CASE DEFAULT
1003  cpabort("Unknown orbital center, try again...")
1004  END SELECT
1005  END DO
1006  !
1007  !
1008  !
1009  ! block the states for each center
1010  ALLOCATE (current_env%psi0_order(nspins))
1011  DO ispin = 1, nspins
1012  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1013  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  context=mo_coeff%matrix_struct%context)
1017  CALL cp_fm_create(current_env%psi0_order(ispin), tmp_fm_struct)
1018  CALL cp_fm_struct_release(tmp_fm_struct)
1019  CALL cp_fm_set_all(current_env%psi0_order(ispin), 0.0_dp)
1020  END DO
1021  !
1022  !
1023  DO ispin = 1, nspins
1024  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
1025  jstate = 0
1026  DO icenter = 1, current_env%nbr_center(ispin)
1027  DO i = current_env%center_list(ispin)%array(1, icenter), &
1028  current_env%center_list(ispin)%array(1, icenter + 1) - 1
1029  !
1030  istate = current_env%center_list(ispin)%array(2, i)
1031  jstate = jstate + 1
1032  !
1033  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  current_env%center_list(ispin)%array(2, i) = jstate
1037  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  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  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  IF (output_unit > 0) THEN
1054  WRITE (output_unit, *) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING '
1055  WRITE (output_unit, *) 'orbital center shifted to match the '
1056  WRITE (output_unit, *) 'center of the box (L/2 L/2 L/2) (superdirty...)'
1057  WRITE (output_unit, *) 'WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING '
1058  END IF
1059  DO ispin = 1, nspins
1060  DO istate = 1, current_env%nbr_center(ispin)
1061  IF (current_env%centers_set(ispin)%array(1, istate) .LE. 0.0_dp) THEN
1062  current_env%centers_set(ispin)%array(1, istate) = &
1063  current_env%centers_set(ispin)%array(1, istate) + cell%hmat(1, 1)
1064  END IF
1065  IF (current_env%centers_set(ispin)%array(2, istate) .LE. 0.0_dp) THEN
1066  current_env%centers_set(ispin)%array(2, istate) = &
1067  current_env%centers_set(ispin)%array(2, istate) + cell%hmat(2, 2)
1068  END IF
1069  IF (current_env%centers_set(ispin)%array(3, istate) .LE. 0.0_dp) THEN
1070  current_env%centers_set(ispin)%array(3, istate) = &
1071  current_env%centers_set(ispin)%array(3, istate) + cell%hmat(3, 3)
1072  END IF
1073  END DO
1074  END DO
1075  ! printing
1076  DO ispin = 1, nspins
1077  IF (output_unit > 0) THEN
1078  WRITE (output_unit, '(/,T2,A,I2)') "WANNIER CENTERS for spin ", ispin
1079  WRITE (output_unit, '(/,T18,A)') "--------------- Shifted Centers --------------- "
1080  DO istate = 1, current_env%nbr_center(ispin)
1081  WRITE (output_unit, '(T5,A,I6,3F12.6)') &
1082  '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  max_nbr_center = maxval(current_env%nbr_center(1:nspins))
1092  current_env%nao = nao
1093  ALLOCATE (current_env%statetrueindex(3, max_nbr_center, nspins))
1094  current_env%statetrueindex(:, :, :) = 0
1095  IF (.true.) THEN
1096  ALLOCATE (state_done(3, max_nbr_center))
1097  DO ispin = 1, nspins
1098  state_done(:, :) = .false.
1099  current_env%statetrueindex(1, 1, ispin) = 1
1100  center(1) = current_env%centers_set(ispin)%array(1, 1)
1101  center(2) = current_env%centers_set(ispin)%array(2, 1)
1102  center(3) = current_env%centers_set(ispin)%array(3, 1)
1103  state_done(1, 1) = .true.
1104  icount = 1
1105  CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
1106  !
1107  DO idir = 1, 3
1108  ini = 1
1109  IF (idir == 1) ini = 2
1110  DO istate = ini, current_env%nbr_center(ispin)
1111  mdist = huge(0.0_dp)
1112  !
1113  DO istate2 = 1, current_env%nbr_center(ispin)
1114  IF (.NOT. state_done(idir, istate2)) THEN
1115  center2(1) = current_env%centers_set(ispin)%array(1, istate2)
1116  center2(2) = current_env%centers_set(ispin)%array(2, istate2)
1117  center2(3) = current_env%centers_set(ispin)%array(3, istate2)
1118  !
1119  rab = pbc(center, center2, cell)
1120  CALL set_vecp(idir, j, k)
1121  dist = sqrt(rab(j)*rab(j) + rab(k)*rab(k))
1122  !
1123  IF (dist .LT. mdist) THEN
1124  mdist = dist
1125  istate_next = istate2
1126  END IF
1127  END IF
1128  END DO ! istate2
1129  !
1130  icount = icount + 1
1131  state_done(idir, istate_next) = .true.
1132  current_env%statetrueindex(idir, icount, ispin) = istate_next
1133  center(1) = current_env%centers_set(ispin)%array(1, istate_next)
1134  center(2) = current_env%centers_set(ispin)%array(2, istate_next)
1135  center(3) = current_env%centers_set(ispin)%array(3, istate_next)
1136  END DO ! istate
1137  icount = 0
1138  END DO ! idir
1139  END DO
1140  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  IF (output_unit > 0) THEN
1151  WRITE (output_unit, "(T2,A,T60,A)") "CURRENT| Gauge used", &
1152  repeat(' ', 20 - len_trim(current_env%gauge_name))//trim(current_env%gauge_name)
1153  WRITE (output_unit, "(T2,A,T79,L1)") "CURRENT| Use old gauge code", current_env%use_old_gauge_atom
1154  WRITE (output_unit, "(T2,A,T79,L1)") "CURRENT| Compute chi for PBC calculation", current_env%chi_pbc
1155  IF (current_env%gauge .EQ. current_gauge_atom) THEN
1156  WRITE (output_unit, "(T2,A,T68,E12.6)") "CURRENT| Radius for building the gauge", &
1157  current_env%gauge_atom_radius
1158  END IF
1159  WRITE (output_unit, "(T2,A,T60,A)") "CURRENT| Orbital center used", &
1160  repeat(' ', 20 - len_trim(current_env%orb_center_name))//trim(current_env%orb_center_name)
1161  IF (current_env%orb_center .EQ. current_orb_center_common) THEN
1162  WRITE (output_unit, "(T2,A,T50,3F10.6)") "CURRENT| Common center", common_center(1:3)
1163  ELSEIF (current_env%orb_center .EQ. current_orb_center_box) THEN
1164  WRITE (output_unit, "(T2,A,T60,3I5)") "CURRENT| Nbr boxes in each direction", nbox(1:3)
1165  END IF
1166  !
1167  DO ispin = 1, nspins
1168  CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
1169  WRITE (output_unit, '(T2,A,I6,A,I6,A,I2)') 'CURRENT| Compute', &
1170  nstate_list(ispin), ' selected response functions out of', &
1171  nmo, ' for spin ', ispin
1172  !
1173  WRITE (output_unit, '(T2,A,I6,A,I2)') 'CURRENT| There is a total of ', &
1174  current_env%nbr_center(ispin), ' (clustered) center(s) for spin ', ispin
1175  END DO
1176  END IF
1177 
1178  IF (btest(cp_print_key_should_output(logger%iter_info, current_section,&
1179  & "PRINT%RESPONSE_FUNCTION_CUBES"), cp_p_file)) THEN
1180 
1181  NULLIFY (bounds, list)
1182  ncubes = 0
1183  CALL section_vals_val_get(current_section,&
1184  & "PRINT%RESPONSE_FUNCTION_CUBES%CUBES_LU_BOUNDS", &
1185  i_vals=bounds)
1186  ncubes = bounds(2) - bounds(1) + 1
1187  IF (ncubes > 0) THEN
1188  ALLOCATE (current_env%list_cubes(ncubes))
1189  DO ir = 1, ncubes
1190  current_env%list_cubes(ir) = bounds(1) + (ir - 1)
1191  END DO
1192  END IF
1193  IF (.NOT. ASSOCIATED(current_env%list_cubes)) THEN
1194  CALL section_vals_val_get(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%CUBES_LIST", &
1195  n_rep_val=n_rep)
1196  ncubes = 0
1197  DO ir = 1, n_rep
1198  NULLIFY (list)
1199  CALL section_vals_val_get(current_section, "PRINT%RESPONSE_FUNCTION_CUBES%CUBES_LIST", &
1200  i_rep_val=ir, i_vals=list)
1201  IF (ASSOCIATED(list)) THEN
1202  CALL reallocate(current_env%list_cubes, 1, ncubes + SIZE(list))
1203  DO ini = 1, SIZE(list)
1204  current_env%list_cubes(ini + ncubes) = list(ini)
1205  END DO
1206  ncubes = ncubes + SIZE(list)
1207  END IF
1208  END DO ! ir
1209  END IF
1210  IF (.NOT. ASSOCIATED(current_env%list_cubes)) THEN
1211  ALLOCATE (current_env%list_cubes(max_states))
1212  DO ir = 1, max_states
1213  current_env%list_cubes(ir) = ir
1214  END DO
1215  END IF
1216  END IF
1217  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  current_env%psi1_p(nspins, 3), current_env%psi1_rxp(nspins, 3))
1236  IF (current_env%full) THEN
1237  ALLOCATE (current_env%psi1_D(nspins, 3))
1238  END IF
1239  DO ispin = 1, nspins
1240  mo_coeff => current_env%psi0_order(ispin)
1241  NULLIFY (tmp_fm_struct)
1242  CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
1243  ncol_global=current_env%nstates(ispin), &
1244  context=mo_coeff%matrix_struct%context)
1245  DO idir = 1, 3
1246  CALL cp_fm_create(current_env%psi1_p(ispin, idir), tmp_fm_struct)
1247  CALL cp_fm_create(current_env%psi1_rxp(ispin, idir), tmp_fm_struct)
1248  CALL cp_fm_create(current_env%p_psi0(ispin, idir), tmp_fm_struct)
1249  CALL cp_fm_create(current_env%rxp_psi0(ispin, idir), tmp_fm_struct)
1250  IF (current_env%full) THEN
1251  CALL cp_fm_create(current_env%psi1_D(ispin, idir), tmp_fm_struct)
1252  END IF
1253  END DO
1254  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  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1259  ALLOCATE (current_env%jrho1_set(3))
1260  DO idir = 1, 3
1261  NULLIFY (rho_r, rho_g)
1262  ALLOCATE (rho_r(nspins), rho_g(nspins))
1263  DO ispin = 1, nspins
1264  CALL auxbas_pw_pool%create_pw(rho_r(ispin))
1265  CALL pw_zero(rho_r(ispin))
1266  CALL auxbas_pw_pool%create_pw(rho_g(ispin))
1267  CALL pw_zero(rho_g(ispin))
1268  END DO
1269  NULLIFY (current_env%jrho1_set(idir)%rho)
1270  ALLOCATE (current_env%jrho1_set(idir)%rho)
1271  CALL qs_rho_create(current_env%jrho1_set(idir)%rho)
1272  CALL qs_rho_set(current_env%jrho1_set(idir)%rho, &
1273  rho_r=rho_r, rho_g=rho_g)
1274  END DO
1275  !
1276  ! Initialize local current density if GAPW calculation
1277  IF (gapw) THEN
1278  CALL init_jrho_atom_set(jrho1_atom_set, atomic_kind_set, nspins)
1279  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  first_sgf(natom), last_sgf(natom))
1284  current_env%basisfun_center = 0.0_dp
1285  CALL get_particle_set(particle_set, qs_kind_set, &
1286  first_sgf=first_sgf, &
1287  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  DO iatom = 1, natom
1291  DO iao = first_sgf(iatom), last_sgf(iatom)
1292  DO idir = 1, 3
1293  current_env%basisfun_center(idir, iao) = particle_set(iatom)%r(idir)
1294  END DO
1295  END DO
1296  END DO
1297 
1298  DEALLOCATE (first_sgf, last_sgf, state_list)
1299 
1300  current_env%simple_done(1:6) = .false.
1301 
1302  ALLOCATE (current_env%full_done(3*nspins, max_nbr_center))
1303  current_env%full_done = .false.
1304  !
1305  ! allocate pointer for the gauge
1306  ALLOCATE (current_env%rs_gauge(3, pw_env%gridlevel_info%ngrid_levels))
1307 
1308  CALL cp_print_key_finished_output(output_unit, logger, lr_section, "PRINT%PROGRAM_RUN_INFO")
1309  CALL timestop(handle)
1310 
1311  END SUBROUTINE current_env_init
1312 
1313 ! **************************************************************************************************
1314 !> \brief ...
1315 !> \param current_env ...
1316 ! **************************************************************************************************
1317  SUBROUTINE current_env_cleanup(current_env)
1318 
1319  TYPE(current_env_type) :: current_env
1320 
1321  INTEGER :: i, idir, ispin
1322 
1323  CALL cp_fm_release(current_env%psi0_order)
1324  !
1325  !psi1_p
1326  CALL cp_fm_release(current_env%psi1_p)
1327  !
1328  !psi1_rxp
1329  CALL cp_fm_release(current_env%psi1_rxp)
1330  !
1331  !psi1_D
1332  CALL cp_fm_release(current_env%psi1_D)
1333  !
1334  !p_psi0
1335  CALL cp_fm_release(current_env%p_psi0)
1336  !
1337  !rxp_psi0
1338  CALL cp_fm_release(current_env%rxp_psi0)
1339  !
1340  IF (ASSOCIATED(current_env%centers_set)) THEN
1341  DO ispin = 1, SIZE(current_env%centers_set, 1)
1342  DEALLOCATE (current_env%centers_set(ispin)%array)
1343  END DO
1344  DEALLOCATE (current_env%centers_set)
1345  END IF
1346 
1347  IF (ASSOCIATED(current_env%center_list)) THEN
1348  DO ispin = 1, SIZE(current_env%center_list, 1)
1349  DEALLOCATE (current_env%center_list(ispin)%array)
1350  END DO
1351  DEALLOCATE (current_env%center_list)
1352  END IF
1353 
1354  IF (ASSOCIATED(current_env%list_cubes)) THEN
1355  DEALLOCATE (current_env%list_cubes)
1356  END IF
1357  ! Current density on the grid
1358  IF (ASSOCIATED(current_env%jrho1_set)) THEN
1359  DO idir = 1, 3
1360  CALL qs_rho_clear(current_env%jrho1_set(idir)%rho)
1361  DEALLOCATE (current_env%jrho1_set(idir)%rho)
1362  END DO
1363  DEALLOCATE (current_env%jrho1_set)
1364  END IF
1365  !
1366  ! Local current density, atom by atom (only gapw)
1367  IF (ASSOCIATED(current_env%jrho1_atom_set)) THEN
1368  CALL deallocate_jrho_atom_set(current_env%jrho1_atom_set)
1369  END IF
1370 
1371  !fullnmr_done
1372  IF (ASSOCIATED(current_env%full_done)) THEN
1373  DEALLOCATE (current_env%full_done)
1374  END IF
1375  IF (ASSOCIATED(current_env%basisfun_center)) THEN
1376  DEALLOCATE (current_env%basisfun_center)
1377  END IF
1378  IF (ASSOCIATED(current_env%statetrueindex)) THEN
1379  DEALLOCATE (current_env%statetrueindex)
1380  END IF
1381  IF (ASSOCIATED(current_env%selected_states_on_atom_list)) THEN
1382  DEALLOCATE (current_env%selected_states_on_atom_list)
1383  END IF
1384  ! give back the gauge
1385  IF (ASSOCIATED(current_env%rs_gauge)) THEN
1386  DO idir = 1, 3
1387  DO i = 1, SIZE(current_env%rs_gauge, 2)
1388  CALL rs_grid_release(current_env%rs_gauge(idir, i))
1389  END DO
1390  END DO
1391  DEALLOCATE (current_env%rs_gauge)
1392  END IF
1393 
1394  ! give back the buf
1395  IF (ASSOCIATED(current_env%rs_buf)) THEN
1396  DO i = 1, SIZE(current_env%rs_buf)
1397  CALL rs_grid_release(current_env%rs_buf(i))
1398  END DO
1399  DEALLOCATE (current_env%rs_buf)
1400  END IF
1401 
1402  END SUBROUTINE current_env_cleanup
1403 
1404 END MODULE qs_linres_current_utils
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1816
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:101
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::current_orb_center_wannier
integer, parameter, public current_orb_center_wannier
Definition: input_constants.F:501
input_constants::current_gauge_atom
integer, parameter, public current_gauge_atom
Definition: input_constants.F:501
input_constants::current_gauge_r
integer, parameter, public current_gauge_r
Definition: input_constants.F:501
input_constants::current_gauge_r_and_step_func
integer, parameter, public current_gauge_r_and_step_func
Definition: input_constants.F:501
input_constants::current_orb_center_box
integer, parameter, public current_orb_center_box
Definition: input_constants.F:501
input_constants::current_orb_center_common
integer, parameter, public current_orb_center_common
Definition: input_constants.F:501
input_constants::current_orb_center_atom
integer, parameter, public current_orb_center_atom
Definition: input_constants.F:501
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition: input_constants.F:513
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition: input_section_types.F:1005
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_path_length
integer, parameter, public default_path_length
Definition: kinds.F:58
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_linres_current_utils
Chemical shift calculation by dfpt Initialization of the nmr_env, creation of the special neighbor li...
Definition: qs_linres_current_utils.F:22
qs_linres_current_utils::current_response
subroutine, public current_response(current_env, p_env, qs_env)
...
Definition: qs_linres_current_utils.F:128
qs_linres_current_utils::current_env_init
subroutine, public current_env_init(current_env, qs_env)
...
Definition: qs_linres_current_utils.F:567
qs_linres_current_utils::current_env_cleanup
subroutine, public current_env_cleanup(current_env)
...
Definition: qs_linres_current_utils.F:1318
qs_linres_methods
localize wavefunctions linear response scf
Definition: qs_linres_methods.F:15
qs_linres_methods::linres_write_restart
subroutine, public linres_write_restart(qs_env, linres_section, vec, ivec, tag, ind)
...
Definition: qs_linres_methods.F:808
qs_linres_methods::linres_read_restart
subroutine, public linres_read_restart(qs_env, linres_section, vec, ivec, tag, ind)
...
Definition: qs_linres_methods.F:922
qs_linres_methods::linres_solver
subroutine, public linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
scf loop to optimize the first order wavefunctions (psi1) given a perturbation as an operator applied...
Definition: qs_linres_methods.F:216
qs_linres_op
Calculate the operators p rxp and D needed in the optimization of the different contribution of the f...
Definition: qs_linres_op.F:18
qs_linres_op::set_vecp
subroutine, public set_vecp(i1, i2, i3)
...
Definition: qs_linres_op.F:1155
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_linres_types::init_jrho_atom_set
subroutine, public init_jrho_atom_set(jrho1_atom_set, atomic_kind_set, nspins)
...
Definition: qs_linres_types.F:1088
qs_linres_types::set_current_env
subroutine, public set_current_env(current_env, jrho1_atom_set, jrho1_set)
...
Definition: qs_linres_types.F:609
qs_linres_types::deallocate_jrho_atom_set
subroutine, public deallocate_jrho_atom_set(jrho_atom_set)
...
Definition: qs_linres_types.F:783
qs_linres_types::get_current_env
subroutine, public get_current_env(current_env, simple_done, simple_converged, full_done, nao, nstates, gauge, list_cubes, statetrueindex, gauge_name, basisfun_center, nbr_center, center_list, centers_set, psi1_p, psi1_rxp, psi1_D, p_psi0, rxp_psi0, jrho1_atom_set, jrho1_set, chi_tensor, chi_tensor_loc, gauge_atom_radius, rs_gauge, use_old_gauge_atom, chi_pbc, psi0_order)
...
Definition: qs_linres_types.F:426
qs_loc_methods
Driver for the localization that should be general for all the methods available and all the definiti...
Definition: qs_loc_methods.F:21
qs_loc_methods::qs_print_cubes
subroutine, public qs_print_cubes(qs_env, mo_coeff, nstates, state_list, centers, print_key, root, ispin, idir, state0, file_position)
write the cube files for a set of selected states
Definition: qs_loc_methods.F:784
qs_loc_types
New version of the module for the localization of the molecular orbitals This should be able to use d...
Definition: qs_loc_types.F:25
qs_loc_types::get_qs_loc_env
subroutine, public get_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
Definition: qs_loc_types.F:317
qs_matrix_pools
wrapper for the pools of matrixes
Definition: qs_matrix_pools.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_operators_ao
Definition: qs_operators_ao.F:13
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition: qs_operators_ao.F:74
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition: qs_p_env_types.F:14
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_set
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
qs_rho_types::qs_rho_create
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
Definition: qs_rho_types.F:99
qs_rho_types::qs_rho_clear
subroutine, public qs_rho_clear(rho_struct)
Deallocates all components, without deallocating rho_struct itself.
Definition: qs_rho_types.F:125
realspace_grid_types
Definition: realspace_grid_types.F:18
realspace_grid_types::rs_grid_release
subroutine, public rs_grid_release(rs_grid)
releases the given rs grid (see doc/ReferenceCounting.html)
Definition: realspace_grid_types.F:2016
scf_control_types
parameters that control an scf iteration
Definition: scf_control_types.F:17