LCOV - code coverage report
Current view: top level - src - qs_linres_current_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.9 % 630 585
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Chemical shift calculation by dfpt
      10              : !>      Initialization of the nmr_env, creation of the special neighbor lists
      11              : !>      Perturbation Hamiltonians by application of the p and rxp oprtators to  psi0
      12              : !>      Write output
      13              : !>      Deallocate everything
      14              : !> \note
      15              : !>      The psi0 should be localized
      16              : !>      the Sebastiani method works within the assumption that the orbitals are
      17              : !>      completely contained in the simulation box
      18              : !> \par History
      19              : !>       created 07-2005 [MI]
      20              : !> \author MI
      21              : ! **************************************************************************************************
      22              : MODULE qs_linres_current_utils
      23              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      24              :    USE cell_types,                      ONLY: cell_type,&
      25              :                                               pbc
      26              :    USE cp_array_utils,                  ONLY: cp_2d_i_p_type,&
      27              :                                               cp_2d_r_p_type
      28              :    USE cp_control_types,                ONLY: dft_control_type
      29              :    USE cp_dbcsr_api,                    ONLY: dbcsr_convert_offsets_to_sizes,&
      30              :                                               dbcsr_copy,&
      31              :                                               dbcsr_create,&
      32              :                                               dbcsr_distribution_type,&
      33              :                                               dbcsr_p_type,&
      34              :                                               dbcsr_set,&
      35              :                                               dbcsr_type_antisymmetric
      36              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      37              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      38              :                                               dbcsr_allocate_matrix_set,&
      39              :                                               dbcsr_deallocate_matrix_set
      40              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      41              :                                               cp_fm_scale_and_add
      42              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      43              :                                               cp_fm_struct_release,&
      44              :                                               cp_fm_struct_type
      45              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      46              :                                               cp_fm_release,&
      47              :                                               cp_fm_set_all,&
      48              :                                               cp_fm_to_fm,&
      49              :                                               cp_fm_type
      50              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      51              :                                               cp_logger_type,&
      52              :                                               cp_to_string
      53              :    USE cp_output_handling,              ONLY: cp_p_file,&
      54              :                                               cp_print_key_finished_output,&
      55              :                                               cp_print_key_should_output,&
      56              :                                               cp_print_key_unit_nr
      57              :    USE input_constants,                 ONLY: current_gauge_atom,&
      58              :                                               current_gauge_r,&
      59              :                                               current_gauge_r_and_step_func,&
      60              :                                               current_orb_center_atom,&
      61              :                                               current_orb_center_box,&
      62              :                                               current_orb_center_common,&
      63              :                                               current_orb_center_wannier,&
      64              :                                               ot_precond_full_all
      65              :    USE input_section_types,             ONLY: section_get_lval,&
      66              :                                               section_vals_get_subs_vals,&
      67              :                                               section_vals_type,&
      68              :                                               section_vals_val_get
      69              :    USE kinds,                           ONLY: default_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 :                            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         6644 :                            mdist = dist
    1125         6644 :                            istate_next = istate2
    1126              :                         END IF
    1127              :                      END IF
    1128              :                   END DO ! istate2
    1129              :                   !
    1130         3500 :                   icount = icount + 1
    1131         3500 :                   state_done(idir, istate_next) = .TRUE.
    1132         3500 :                   current_env%statetrueindex(idir, icount, ispin) = istate_next
    1133         3500 :                   center(1) = current_env%centers_set(ispin)%array(1, istate_next)
    1134         3500 :                   center(2) = current_env%centers_set(ispin)%array(2, istate_next)
    1135         4250 :                   center(3) = current_env%centers_set(ispin)%array(3, istate_next)
    1136              :                END DO ! istate
    1137         1000 :                icount = 0
    1138              :             END DO ! idir
    1139              :          END DO
    1140          174 :          DEALLOCATE (state_done)
    1141              :       ELSE
    1142              :          DO ispin = 1, nspins
    1143              :             DO icenter = 1, current_env%nbr_center(ispin)
    1144              :                current_env%statetrueindex(:, icenter, ispin) = icenter
    1145              :             END DO
    1146              :          END DO
    1147              :       END IF
    1148              :       !
    1149              :       !
    1150          174 :       IF (output_unit > 0) THEN
    1151           87 :          WRITE (output_unit, "(T2,A,T60,A)") "CURRENT| Gauge used", &
    1152          507 :             REPEAT(' ', 20 - LEN_TRIM(current_env%gauge_name))//TRIM(current_env%gauge_name)
    1153           87 :          WRITE (output_unit, "(T2,A,T79,L1)") "CURRENT| Use old gauge code", current_env%use_old_gauge_atom
    1154           87 :          WRITE (output_unit, "(T2,A,T79,L1)") "CURRENT| Compute chi for PBC calculation", current_env%chi_pbc
    1155           87 :          IF (current_env%gauge .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          696 :       ALLOCATE (current_env%jrho1_set(3))
    1260          696 :       DO idir = 1, 3
    1261          522 :          NULLIFY (rho_r, rho_g)
    1262         4110 :          ALLOCATE (rho_r(nspins), rho_g(nspins))
    1263         1272 :          DO ispin = 1, nspins
    1264          750 :             CALL auxbas_pw_pool%create_pw(rho_r(ispin))
    1265          750 :             CALL pw_zero(rho_r(ispin))
    1266          750 :             CALL auxbas_pw_pool%create_pw(rho_g(ispin))
    1267         1272 :             CALL pw_zero(rho_g(ispin))
    1268              :          END DO
    1269          522 :          NULLIFY (current_env%jrho1_set(idir)%rho)
    1270          522 :          ALLOCATE (current_env%jrho1_set(idir)%rho)
    1271          522 :          CALL qs_rho_create(current_env%jrho1_set(idir)%rho)
    1272              :          CALL qs_rho_set(current_env%jrho1_set(idir)%rho, &
    1273          696 :                          rho_r=rho_r, rho_g=rho_g)
    1274              :       END DO
    1275              :       !
    1276              :       ! Initialize local current density if GAPW calculation
    1277          174 :       IF (gapw) THEN
    1278           96 :          CALL init_jrho_atom_set(jrho1_atom_set, atomic_kind_set, nspins)
    1279           96 :          CALL set_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set)
    1280              :       END IF
    1281              :       !
    1282              :       ALLOCATE (current_env%basisfun_center(3, current_env%nao), &
    1283         1044 :                 first_sgf(natom), last_sgf(natom))
    1284        15886 :       current_env%basisfun_center = 0.0_dp
    1285              :       CALL get_particle_set(particle_set, qs_kind_set, &
    1286              :                             first_sgf=first_sgf, &
    1287          174 :                             last_sgf=last_sgf)
    1288              :       !Build 3 arrays where for each contracted basis function
    1289              :       !the x y and z coordinates of the center are given
    1290          710 :       DO iatom = 1, natom
    1291         4638 :          DO iao = first_sgf(iatom), last_sgf(iatom)
    1292        16248 :             DO idir = 1, 3
    1293        15712 :                current_env%basisfun_center(idir, iao) = particle_set(iatom)%r(idir)
    1294              :             END DO
    1295              :          END DO
    1296              :       END DO
    1297              : 
    1298          174 :       DEALLOCATE (first_sgf, last_sgf, state_list)
    1299              : 
    1300         1218 :       current_env%simple_done(1:6) = .FALSE.
    1301              : 
    1302          696 :       ALLOCATE (current_env%full_done(3*nspins, max_nbr_center))
    1303         5052 :       current_env%full_done = .FALSE.
    1304              :       !
    1305              :       ! allocate pointer for the gauge
    1306         5900 :       ALLOCATE (current_env%rs_gauge(3, pw_env%gridlevel_info%ngrid_levels))
    1307              : 
    1308          174 :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, "PRINT%PROGRAM_RUN_INFO")
    1309          174 :       CALL timestop(handle)
    1310              : 
    1311          870 :    END SUBROUTINE current_env_init
    1312              : 
    1313              : ! **************************************************************************************************
    1314              : !> \brief ...
    1315              : !> \param current_env ...
    1316              : ! **************************************************************************************************
    1317          348 :    SUBROUTINE current_env_cleanup(current_env)
    1318              : 
    1319              :       TYPE(current_env_type)                             :: current_env
    1320              : 
    1321              :       INTEGER                                            :: i, idir, ispin
    1322              : 
    1323          348 :       CALL cp_fm_release(current_env%psi0_order)
    1324              :       !
    1325              :       !psi1_p
    1326          348 :       CALL cp_fm_release(current_env%psi1_p)
    1327              :       !
    1328              :       !psi1_rxp
    1329          348 :       CALL cp_fm_release(current_env%psi1_rxp)
    1330              :       !
    1331              :       !psi1_D
    1332          348 :       CALL cp_fm_release(current_env%psi1_D)
    1333              :       !
    1334              :       !p_psi0
    1335          348 :       CALL cp_fm_release(current_env%p_psi0)
    1336              :       !
    1337              :       !rxp_psi0
    1338          348 :       CALL cp_fm_release(current_env%rxp_psi0)
    1339              :       !
    1340          348 :       IF (ASSOCIATED(current_env%centers_set)) THEN
    1341          424 :          DO ispin = 1, SIZE(current_env%centers_set, 1)
    1342          424 :             DEALLOCATE (current_env%centers_set(ispin)%array)
    1343              :          END DO
    1344          174 :          DEALLOCATE (current_env%centers_set)
    1345              :       END IF
    1346              : 
    1347          348 :       IF (ASSOCIATED(current_env%center_list)) THEN
    1348          424 :          DO ispin = 1, SIZE(current_env%center_list, 1)
    1349          424 :             DEALLOCATE (current_env%center_list(ispin)%array)
    1350              :          END DO
    1351          174 :          DEALLOCATE (current_env%center_list)
    1352              :       END IF
    1353              : 
    1354          348 :       IF (ASSOCIATED(current_env%list_cubes)) THEN
    1355            0 :          DEALLOCATE (current_env%list_cubes)
    1356              :       END IF
    1357              :       ! Current density on the grid
    1358          348 :       IF (ASSOCIATED(current_env%jrho1_set)) THEN
    1359          696 :          DO idir = 1, 3
    1360          522 :             CALL qs_rho_clear(current_env%jrho1_set(idir)%rho)
    1361          696 :             DEALLOCATE (current_env%jrho1_set(idir)%rho)
    1362              :          END DO
    1363          174 :          DEALLOCATE (current_env%jrho1_set)
    1364              :       END IF
    1365              :       !
    1366              :       ! Local current density, atom by atom (only gapw)
    1367          348 :       IF (ASSOCIATED(current_env%jrho1_atom_set)) THEN
    1368           96 :          CALL deallocate_jrho_atom_set(current_env%jrho1_atom_set)
    1369              :       END IF
    1370              : 
    1371              :       !fullnmr_done
    1372          348 :       IF (ASSOCIATED(current_env%full_done)) THEN
    1373          174 :          DEALLOCATE (current_env%full_done)
    1374              :       END IF
    1375          348 :       IF (ASSOCIATED(current_env%basisfun_center)) THEN
    1376          174 :          DEALLOCATE (current_env%basisfun_center)
    1377              :       END IF
    1378          348 :       IF (ASSOCIATED(current_env%statetrueindex)) THEN
    1379          174 :          DEALLOCATE (current_env%statetrueindex)
    1380              :       END IF
    1381          348 :       IF (ASSOCIATED(current_env%selected_states_on_atom_list)) THEN
    1382            8 :          DEALLOCATE (current_env%selected_states_on_atom_list)
    1383              :       END IF
    1384              :       ! give back the gauge
    1385          348 :       IF (ASSOCIATED(current_env%rs_gauge)) THEN
    1386          696 :          DO idir = 1, 3
    1387         2772 :             DO i = 1, SIZE(current_env%rs_gauge, 2)
    1388         2598 :                CALL rs_grid_release(current_env%rs_gauge(idir, i))
    1389              :             END DO
    1390              :          END DO
    1391          174 :          DEALLOCATE (current_env%rs_gauge)
    1392              :       END IF
    1393              : 
    1394              :       ! give back the buf
    1395          348 :       IF (ASSOCIATED(current_env%rs_buf)) THEN
    1396          136 :          DO i = 1, SIZE(current_env%rs_buf)
    1397          136 :             CALL rs_grid_release(current_env%rs_buf(i))
    1398              :          END DO
    1399           28 :          DEALLOCATE (current_env%rs_buf)
    1400              :       END IF
    1401              : 
    1402          348 :    END SUBROUTINE current_env_cleanup
    1403              : 
    1404              : END MODULE qs_linres_current_utils
        

Generated by: LCOV version 2.0-1