LCOV - code coverage report
Current view: top level - src - qs_rho_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 445 633 70.3 %
Date: 2024-03-28 07:31:50 Functions: 8 8 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief methods of the rho structure (defined in qs_rho_types)
      10             : !> \par History
      11             : !>      08.2002 created [fawzi]
      12             : !>      08.2014 kpoints [JGH]
      13             : !> \author Fawzi Mohamed
      14             : ! **************************************************************************************************
      15             : MODULE qs_rho_methods
      16             :    USE admm_types,                      ONLY: get_admm_env
      17             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      18             :    USE cp_control_types,                ONLY: dft_control_type
      19             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      20             :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      21             :                                               dbcsr_deallocate_matrix_set
      22             :    USE cp_log_handling,                 ONLY: cp_to_string
      23             :    USE dbcsr_api,                       ONLY: &
      24             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_p_type, dbcsr_scale, dbcsr_set, dbcsr_type, &
      25             :         dbcsr_type_antisymmetric, dbcsr_type_symmetric
      26             :    USE kinds,                           ONLY: default_string_length,&
      27             :                                               dp
      28             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      29             :                                               kpoint_type
      30             :    USE lri_environment_methods,         ONLY: calculate_lri_densities
      31             :    USE lri_environment_types,           ONLY: lri_density_type,&
      32             :                                               lri_environment_type
      33             :    USE message_passing,                 ONLY: mp_para_env_type
      34             :    USE pw_env_types,                    ONLY: pw_env_get,&
      35             :                                               pw_env_type
      36             :    USE pw_methods,                      ONLY: pw_axpy,&
      37             :                                               pw_copy,&
      38             :                                               pw_scale,&
      39             :                                               pw_zero
      40             :    USE pw_pool_types,                   ONLY: pw_pool_type
      41             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      42             :                                               pw_r3d_rs_type
      43             :    USE qs_collocate_density,            ONLY: calculate_drho_elec,&
      44             :                                               calculate_rho_elec
      45             :    USE qs_environment_types,            ONLY: get_qs_env,&
      46             :                                               qs_environment_type,&
      47             :                                               set_qs_env
      48             :    USE qs_kind_types,                   ONLY: qs_kind_type
      49             :    USE qs_ks_types,                     ONLY: get_ks_env,&
      50             :                                               qs_ks_env_type
      51             :    USE qs_local_rho_types,              ONLY: local_rho_type
      52             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      53             :    USE qs_oce_types,                    ONLY: oce_matrix_type
      54             :    USE qs_rho_atom_methods,             ONLY: calculate_rho_atom_coeff
      55             :    USE qs_rho_atom_types,               ONLY: rho_atom_type
      56             :    USE qs_rho_types,                    ONLY: qs_rho_clear,&
      57             :                                               qs_rho_get,&
      58             :                                               qs_rho_set,&
      59             :                                               qs_rho_type
      60             :    USE ri_environment_methods,          ONLY: calculate_ri_densities
      61             :    USE task_list_types,                 ONLY: task_list_type
      62             : #include "./base/base_uses.f90"
      63             : 
      64             :    IMPLICIT NONE
      65             :    PRIVATE
      66             : 
      67             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      68             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_methods'
      69             : 
      70             :    PUBLIC :: qs_rho_update_rho, qs_rho_update_tddfpt, &
      71             :              qs_rho_rebuild, qs_rho_copy, qs_rho_scale_and_add
      72             :    PUBLIC :: duplicate_rho_type, allocate_rho_ao_imag_from_real
      73             : 
      74             : CONTAINS
      75             : 
      76             : ! **************************************************************************************************
      77             : !> \brief rebuilds rho (if necessary allocating and initializing it)
      78             : !> \param rho the rho type to rebuild (defaults to qs_env%rho)
      79             : !> \param qs_env the environment to which rho belongs
      80             : !> \param rebuild_ao if it is necessary to rebuild rho_ao. Defaults to true.
      81             : !> \param rebuild_grids if it in necessary to rebuild rho_r and rho_g.
      82             : !>        Defaults to false.
      83             : !> \param admm (use aux_fit basis)
      84             : !> \param pw_env_external external plane wave environment
      85             : !> \par History
      86             : !>      11.2002 created replacing qs_rho_create and qs_env_rebuild_rho[fawzi]
      87             : !> \author Fawzi Mohamed
      88             : !> \note
      89             : !>      needs updated  pw pools, s, s_mstruct and h in qs_env.
      90             : !>      The use of p to keep the structure of h (needed for the forces)
      91             : !>      is ugly and should be removed.
      92             : !>      Change so that it does not allocate a subcomponent if it is not
      93             : !>      associated and not requested?
      94             : ! **************************************************************************************************
      95       49108 :    SUBROUTINE qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
      96             :       TYPE(qs_rho_type), INTENT(INOUT)                   :: rho
      97             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      98             :       LOGICAL, INTENT(in), OPTIONAL                      :: rebuild_ao, rebuild_grids, admm
      99             :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
     100             : 
     101             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'qs_rho_rebuild'
     102             : 
     103             :       CHARACTER(LEN=default_string_length)               :: headline
     104             :       INTEGER                                            :: handle, i, ic, j, nimg, nspins
     105             :       LOGICAL                                            :: do_kpoints, my_admm, my_rebuild_ao, &
     106             :                                                             my_rebuild_grids, rho_ao_is_complex
     107       24554 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_r
     108       24554 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_kp, rho_ao_im_kp, rho_ao_kp
     109             :       TYPE(dbcsr_type), POINTER                          :: refmatrix, tmatrix
     110             :       TYPE(dft_control_type), POINTER                    :: dft_control
     111             :       TYPE(kpoint_type), POINTER                         :: kpoints
     112             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     113       24554 :          POINTER                                         :: sab_orb
     114       24554 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g, tau_g
     115       24554 :       TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER     :: drho_g
     116             :       TYPE(pw_env_type), POINTER                         :: pw_env
     117             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     118       24554 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r, tau_r
     119       24554 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER     :: drho_r
     120             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_r_sccs
     121             : 
     122       24554 :       CALL timeset(routineN, handle)
     123             : 
     124       24554 :       NULLIFY (pw_env, auxbas_pw_pool, matrix_s_kp, dft_control)
     125       24554 :       NULLIFY (tot_rho_r, rho_ao_kp, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g, rho_ao_im_kp)
     126       24554 :       NULLIFY (rho_r_sccs)
     127       24554 :       NULLIFY (sab_orb)
     128       24554 :       my_rebuild_ao = .TRUE.
     129       24554 :       my_rebuild_grids = .TRUE.
     130       24554 :       my_admm = .FALSE.
     131       24554 :       IF (PRESENT(rebuild_ao)) my_rebuild_ao = rebuild_ao
     132       24554 :       IF (PRESENT(rebuild_grids)) my_rebuild_grids = rebuild_grids
     133       24554 :       IF (PRESENT(admm)) my_admm = admm
     134             : 
     135             :       CALL get_qs_env(qs_env, &
     136             :                       kpoints=kpoints, &
     137             :                       do_kpoints=do_kpoints, &
     138             :                       pw_env=pw_env, &
     139       24554 :                       dft_control=dft_control)
     140       24554 :       IF (PRESENT(pw_env_external)) &
     141         718 :          pw_env => pw_env_external
     142             : 
     143       24554 :       nimg = dft_control%nimages
     144             : 
     145       24554 :       IF (my_admm) THEN
     146        1790 :          CALL get_admm_env(qs_env%admm_env, sab_aux_fit=sab_orb, matrix_s_aux_fit_kp=matrix_s_kp)
     147             :       ELSE
     148       22764 :          CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp)
     149             : 
     150       22764 :          IF (do_kpoints) THEN
     151         766 :             CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
     152             :          ELSE
     153       21998 :             CALL get_qs_env(qs_env, sab_orb=sab_orb)
     154             :          END IF
     155             :       END IF
     156       24554 :       refmatrix => matrix_s_kp(1, 1)%matrix
     157             : 
     158       24554 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     159       24554 :       nspins = dft_control%nspins
     160             : 
     161             :       CALL qs_rho_get(rho, &
     162             :                       tot_rho_r=tot_rho_r, &
     163             :                       rho_ao_kp=rho_ao_kp, &
     164             :                       rho_ao_im_kp=rho_ao_im_kp, &
     165             :                       rho_r=rho_r, &
     166             :                       rho_g=rho_g, &
     167             :                       drho_r=drho_r, &
     168             :                       drho_g=drho_g, &
     169             :                       tau_r=tau_r, &
     170             :                       tau_g=tau_g, &
     171             :                       rho_r_sccs=rho_r_sccs, &
     172       24554 :                       complex_rho_ao=rho_ao_is_complex)
     173             : 
     174       24554 :       IF (.NOT. ASSOCIATED(tot_rho_r)) THEN
     175       32028 :          ALLOCATE (tot_rho_r(nspins))
     176       23439 :          tot_rho_r = 0.0_dp
     177       10676 :          CALL qs_rho_set(rho, tot_rho_r=tot_rho_r)
     178             :       END IF
     179             : 
     180             :       ! rho_ao
     181       24554 :       IF (my_rebuild_ao .OR. (.NOT. ASSOCIATED(rho_ao_kp))) THEN
     182       23104 :          IF (ASSOCIATED(rho_ao_kp)) &
     183       13868 :             CALL dbcsr_deallocate_matrix_set(rho_ao_kp)
     184             :          ! Create a new density matrix set
     185       23104 :          CALL dbcsr_allocate_matrix_set(rho_ao_kp, nspins, nimg)
     186       23104 :          CALL qs_rho_set(rho, rho_ao_kp=rho_ao_kp)
     187       49338 :          DO i = 1, nspins
     188      149982 :             DO ic = 1, nimg
     189      100644 :                IF (nspins > 1) THEN
     190       26540 :                   IF (i == 1) THEN
     191       13270 :                      headline = "DENSITY MATRIX FOR ALPHA SPIN"
     192             :                   ELSE
     193       13270 :                      headline = "DENSITY MATRIX FOR BETA SPIN"
     194             :                   END IF
     195             :                ELSE
     196       74104 :                   headline = "DENSITY MATRIX"
     197             :                END IF
     198      100644 :                ALLOCATE (rho_ao_kp(i, ic)%matrix)
     199      100644 :                tmatrix => rho_ao_kp(i, ic)%matrix
     200             :                CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=TRIM(headline), &
     201      100644 :                                  matrix_type=dbcsr_type_symmetric, nze=0)
     202      100644 :                CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
     203      126878 :                CALL dbcsr_set(tmatrix, 0.0_dp)
     204             :             END DO
     205             :          END DO
     206       23104 :          IF (rho_ao_is_complex) THEN
     207         328 :             IF (ASSOCIATED(rho_ao_im_kp)) THEN
     208         328 :                CALL dbcsr_deallocate_matrix_set(rho_ao_im_kp)
     209             :             END IF
     210         328 :             CALL dbcsr_allocate_matrix_set(rho_ao_im_kp, nspins, nimg)
     211         328 :             CALL qs_rho_set(rho, rho_ao_im_kp=rho_ao_im_kp)
     212         720 :             DO i = 1, nspins
     213        1112 :                DO ic = 1, nimg
     214         392 :                   IF (nspins > 1) THEN
     215         128 :                      IF (i == 1) THEN
     216          64 :                         headline = "IMAGINARY PART OF DENSITY MATRIX FOR ALPHA SPIN"
     217             :                      ELSE
     218          64 :                         headline = "IMAGINARY PART OF DENSITY MATRIX FOR BETA SPIN"
     219             :                      END IF
     220             :                   ELSE
     221         264 :                      headline = "IMAGINARY PART OF DENSITY MATRIX"
     222             :                   END IF
     223         392 :                   ALLOCATE (rho_ao_im_kp(i, ic)%matrix)
     224         392 :                   tmatrix => rho_ao_im_kp(i, ic)%matrix
     225             :                   CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=TRIM(headline), &
     226         392 :                                     matrix_type=dbcsr_type_antisymmetric, nze=0)
     227         392 :                   CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb)
     228         784 :                   CALL dbcsr_set(tmatrix, 0.0_dp)
     229             :                END DO
     230             :             END DO
     231             :          END IF
     232             :       END IF
     233             : 
     234             :       ! rho_r
     235       24554 :       IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_r)) THEN
     236       24554 :          IF (ASSOCIATED(rho_r)) THEN
     237       28899 :             DO i = 1, SIZE(rho_r)
     238       28899 :                CALL rho_r(i)%release()
     239             :             END DO
     240       13878 :             DEALLOCATE (rho_r)
     241             :          END IF
     242      101446 :          ALLOCATE (rho_r(nspins))
     243       24554 :          CALL qs_rho_set(rho, rho_r=rho_r)
     244       52338 :          DO i = 1, nspins
     245       52338 :             CALL auxbas_pw_pool%create_pw(rho_r(i))
     246             :          END DO
     247             :       END IF
     248             : 
     249             :       ! rho_g
     250       24554 :       IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_g)) THEN
     251       24554 :          IF (ASSOCIATED(rho_g)) THEN
     252       28899 :             DO i = 1, SIZE(rho_g)
     253       28899 :                CALL rho_g(i)%release()
     254             :             END DO
     255       13878 :             DEALLOCATE (rho_g)
     256             :          END IF
     257      101446 :          ALLOCATE (rho_g(nspins))
     258       24554 :          CALL qs_rho_set(rho, rho_g=rho_g)
     259       52338 :          DO i = 1, nspins
     260       52338 :             CALL auxbas_pw_pool%create_pw(rho_g(i))
     261             :          END DO
     262             :       END IF
     263             : 
     264             :       ! SCCS
     265       24554 :       IF (dft_control%do_sccs) THEN
     266          12 :          IF (my_rebuild_grids .OR. (.NOT. ASSOCIATED(rho_r_sccs))) THEN
     267          12 :             IF (ASSOCIATED(rho_r_sccs)) THEN
     268           2 :                CALL rho_r_sccs%release()
     269           2 :                DEALLOCATE (rho_r_sccs)
     270             :             END IF
     271          12 :             ALLOCATE (rho_r_sccs)
     272          12 :             CALL qs_rho_set(rho, rho_r_sccs=rho_r_sccs)
     273          12 :             CALL auxbas_pw_pool%create_pw(rho_r_sccs)
     274          12 :             CALL pw_zero(rho_r_sccs)
     275             :          END IF
     276             :       END IF
     277             : 
     278             :       ! allocate drho_r and drho_g if xc_deriv_collocate
     279       24554 :       IF (dft_control%drho_by_collocation) THEN
     280             :          ! drho_r
     281           0 :          IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_r)) THEN
     282           0 :             IF (ASSOCIATED(drho_r)) THEN
     283           0 :                DO j = 1, SIZE(drho_r, 2)
     284           0 :                   DO i = 1, SIZE(drho_r, 1)
     285           0 :                      CALL drho_r(i, j)%release()
     286             :                   END DO
     287             :                END DO
     288           0 :                DEALLOCATE (drho_r)
     289             :             END IF
     290           0 :             ALLOCATE (drho_r(3, nspins))
     291           0 :             CALL qs_rho_set(rho, drho_r=drho_r)
     292           0 :             DO j = 1, nspins
     293           0 :                DO i = 1, 3
     294           0 :                   CALL auxbas_pw_pool%create_pw(drho_r(i, j))
     295             :                END DO
     296             :             END DO
     297             :          END IF
     298             :          ! drho_g
     299           0 :          IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_g)) THEN
     300           0 :             IF (ASSOCIATED(drho_g)) THEN
     301           0 :                DO j = 1, SIZE(drho_g, 2)
     302           0 :                   DO i = 1, SIZE(drho_r, 1)
     303           0 :                      CALL drho_g(i, j)%release()
     304             :                   END DO
     305             :                END DO
     306           0 :                DEALLOCATE (drho_g)
     307             :             END IF
     308           0 :             ALLOCATE (drho_g(3, nspins))
     309           0 :             CALL qs_rho_set(rho, drho_g=drho_g)
     310           0 :             DO j = 1, nspins
     311           0 :                DO i = 1, 3
     312           0 :                   CALL auxbas_pw_pool%create_pw(drho_g(i, j))
     313             :                END DO
     314             :             END DO
     315             :          END IF
     316             :       END IF
     317             : 
     318             :       ! allocate tau_r and tau_g if use_kinetic_energy_density
     319       24554 :       IF (dft_control%use_kinetic_energy_density) THEN
     320             :          ! tau_r
     321         300 :          IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_r)) THEN
     322         300 :             IF (ASSOCIATED(tau_r)) THEN
     323         118 :                DO i = 1, SIZE(tau_r)
     324         118 :                   CALL tau_r(i)%release()
     325             :                END DO
     326          56 :                DEALLOCATE (tau_r)
     327             :             END IF
     328        1246 :             ALLOCATE (tau_r(nspins))
     329         300 :             CALL qs_rho_set(rho, tau_r=tau_r)
     330         646 :             DO i = 1, nspins
     331         646 :                CALL auxbas_pw_pool%create_pw(tau_r(i))
     332             :             END DO
     333             :          END IF
     334             : 
     335             :          ! tau_g
     336         300 :          IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_g)) THEN
     337         300 :             IF (ASSOCIATED(tau_g)) THEN
     338         118 :                DO i = 1, SIZE(tau_g)
     339         118 :                   CALL tau_g(i)%release()
     340             :                END DO
     341          56 :                DEALLOCATE (tau_g)
     342             :             END IF
     343        1246 :             ALLOCATE (tau_g(nspins))
     344         300 :             CALL qs_rho_set(rho, tau_g=tau_g)
     345         646 :             DO i = 1, nspins
     346         646 :                CALL auxbas_pw_pool%create_pw(tau_g(i))
     347             :             END DO
     348             :          END IF
     349             :       END IF ! use_kinetic_energy_density
     350             : 
     351       24554 :       CALL timestop(handle)
     352             : 
     353       24554 :    END SUBROUTINE qs_rho_rebuild
     354             : 
     355             : ! **************************************************************************************************
     356             : !> \brief updates rho_r and rho_g to the rho%rho_ao.
     357             : !>      if use_kinetic_energy_density also computes tau_r and tau_g
     358             : !>      this works for all ground state and ground state response methods
     359             : !> \param rho_struct the rho structure that should be updated
     360             : !> \param qs_env the qs_env rho_struct refers to
     361             : !>        the integrated charge in r space
     362             : !> \param rho_xc_external ...
     363             : !> \param local_rho_set ...
     364             : !> \param task_list_external external task list
     365             : !> \param task_list_external_soft external task list (soft_version)
     366             : !> \param pw_env_external    external plane wave environment
     367             : !> \param para_env_external  external MPI environment
     368             : !> \par History
     369             : !>      08.2002 created [fawzi]
     370             : !> \author Fawzi Mohamed
     371             : ! **************************************************************************************************
     372      206075 :    SUBROUTINE qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, &
     373             :                                 task_list_external, task_list_external_soft, &
     374             :                                 pw_env_external, para_env_external)
     375             :       TYPE(qs_rho_type), INTENT(INOUT)                   :: rho_struct
     376             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     377             :       TYPE(qs_rho_type), OPTIONAL, POINTER               :: rho_xc_external
     378             :       TYPE(local_rho_type), OPTIONAL, POINTER            :: local_rho_set
     379             :       TYPE(task_list_type), OPTIONAL, POINTER            :: task_list_external, &
     380             :                                                             task_list_external_soft
     381             :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
     382             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env_external
     383             : 
     384      206075 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     385      206075 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     386      206075 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
     387      206075 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
     388             :       TYPE(dft_control_type), POINTER                    :: dft_control
     389             :       TYPE(kpoint_type), POINTER                         :: kpoints
     390             :       TYPE(lri_density_type), POINTER                    :: lri_density
     391             :       TYPE(lri_environment_type), POINTER                :: lri_env
     392             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     393             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     394             : 
     395             :       CALL get_qs_env(qs_env, dft_control=dft_control, &
     396             :                       atomic_kind_set=atomic_kind_set, &
     397      206075 :                       para_env=para_env)
     398      206075 :       IF (PRESENT(para_env_external)) para_env => para_env_external
     399             : 
     400             :       IF (dft_control%qs_control%semi_empirical .OR. &
     401      206075 :           dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
     402             : 
     403       78418 :          CALL qs_rho_set(rho_struct, rho_r_valid=.FALSE., rho_g_valid=.FALSE.)
     404             : 
     405      127657 :       ELSEIF (dft_control%qs_control%lrigpw) THEN
     406         618 :          CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
     407         618 :          CPASSERT(.NOT. dft_control%drho_by_collocation)
     408         618 :          CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     409         618 :          CALL get_qs_env(qs_env, ks_env=ks_env)
     410         618 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
     411         618 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     412         618 :          CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
     413             :          CALL calculate_lri_densities(lri_env, lri_density, qs_env, rho_ao_kp, cell_to_index, &
     414             :                                       lri_rho_struct=rho_struct, &
     415             :                                       atomic_kind_set=atomic_kind_set, &
     416             :                                       para_env=para_env, &
     417         618 :                                       response_density=.FALSE.)
     418         618 :          CALL set_qs_env(qs_env, lri_density=lri_density)
     419         618 :          CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     420             : 
     421      127039 :       ELSEIF (dft_control%qs_control%rigpw) THEN
     422           0 :          CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
     423           0 :          CPASSERT(.NOT. dft_control%drho_by_collocation)
     424           0 :          CALL get_qs_env(qs_env, lri_env=lri_env)
     425           0 :          CALL qs_rho_get(rho_struct, rho_ao=rho_ao)
     426             :          CALL calculate_ri_densities(lri_env, qs_env, rho_ao, &
     427             :                                      lri_rho_struct=rho_struct, &
     428             :                                      atomic_kind_set=atomic_kind_set, &
     429           0 :                                      para_env=para_env)
     430           0 :          CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     431             : 
     432             :       ELSE
     433             :          CALL qs_rho_update_rho_low(rho_struct=rho_struct, qs_env=qs_env, &
     434             :                                     rho_xc_external=rho_xc_external, &
     435             :                                     local_rho_set=local_rho_set, &
     436             :                                     task_list_external=task_list_external, &
     437             :                                     task_list_external_soft=task_list_external_soft, &
     438             :                                     pw_env_external=pw_env_external, &
     439      127039 :                                     para_env_external=para_env_external)
     440             : 
     441             :       END IF
     442             : 
     443      206075 :    END SUBROUTINE qs_rho_update_rho
     444             : 
     445             : ! **************************************************************************************************
     446             : !> \brief updates rho_r and rho_g to the rho%rho_ao.
     447             : !>      if use_kinetic_energy_density also computes tau_r and tau_g
     448             : !> \param rho_struct the rho structure that should be updated
     449             : !> \param qs_env the qs_env rho_struct refers to
     450             : !>        the integrated charge in r space
     451             : !> \param rho_xc_external rho structure for GAPW_XC
     452             : !> \param local_rho_set ...
     453             : !> \param pw_env_external    external plane wave environment
     454             : !> \param task_list_external external task list (use for default and GAPW)
     455             : !> \param task_list_external_soft external task list (soft density for GAPW_XC)
     456             : !> \param para_env_external ...
     457             : !> \par History
     458             : !>      08.2002 created [fawzi]
     459             : !> \author Fawzi Mohamed
     460             : ! **************************************************************************************************
     461      127039 :    SUBROUTINE qs_rho_update_rho_low(rho_struct, qs_env, rho_xc_external, &
     462             :                                     local_rho_set, pw_env_external, &
     463             :                                     task_list_external, task_list_external_soft, &
     464             :                                     para_env_external)
     465             :       TYPE(qs_rho_type), INTENT(INOUT)                   :: rho_struct
     466             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     467             :       TYPE(qs_rho_type), OPTIONAL, POINTER               :: rho_xc_external
     468             :       TYPE(local_rho_type), OPTIONAL, POINTER            :: local_rho_set
     469             :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
     470             :       TYPE(task_list_type), OPTIONAL, POINTER            :: task_list_external, &
     471             :                                                             task_list_external_soft
     472             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env_external
     473             : 
     474             :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_update_rho_low'
     475             : 
     476             :       INTEGER                                            :: handle, img, ispin, nimg, nspins
     477             :       LOGICAL                                            :: gapw, gapw_xc
     478             :       REAL(KIND=dp)                                      :: dum
     479      127039 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_r, tot_rho_r_xc
     480      127039 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     481      127039 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
     482      127039 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp, rho_xc_ao
     483             :       TYPE(dft_control_type), POINTER                    :: dft_control
     484             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     485             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     486      127039 :          POINTER                                         :: sab
     487             :       TYPE(oce_matrix_type), POINTER                     :: oce
     488      127039 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g, rho_xc_g, tau_g, tau_xc_g
     489      127039 :       TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER     :: drho_g, drho_xc_g
     490             :       TYPE(pw_env_type), POINTER                         :: pw_env
     491      127039 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r, rho_xc_r, tau_r, tau_xc_r
     492      127039 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER     :: drho_r, drho_xc_r
     493      127039 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     494             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     495             :       TYPE(qs_rho_type), POINTER                         :: rho_xc
     496      127039 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
     497             :       TYPE(task_list_type), POINTER                      :: task_list
     498             : 
     499      127039 :       CALL timeset(routineN, handle)
     500             : 
     501      127039 :       NULLIFY (dft_control, rho_xc, ks_env, rho_ao, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g)
     502      127039 :       NULLIFY (rho_xc_ao, rho_xc_g, rho_xc_r, drho_xc_g, tau_xc_r, tau_xc_g, tot_rho_r, tot_rho_r_xc)
     503      127039 :       NULLIFY (para_env, pw_env, atomic_kind_set)
     504             : 
     505             :       CALL get_qs_env(qs_env, &
     506             :                       ks_env=ks_env, &
     507             :                       dft_control=dft_control, &
     508      127039 :                       atomic_kind_set=atomic_kind_set)
     509             : 
     510             :       CALL qs_rho_get(rho_struct, &
     511             :                       rho_r=rho_r, &
     512             :                       rho_g=rho_g, &
     513             :                       tot_rho_r=tot_rho_r, &
     514             :                       drho_r=drho_r, &
     515             :                       drho_g=drho_g, &
     516             :                       tau_r=tau_r, &
     517      127039 :                       tau_g=tau_g)
     518             : 
     519             :       CALL get_qs_env(qs_env, task_list=task_list, &
     520      127039 :                       para_env=para_env, pw_env=pw_env)
     521      127039 :       IF (PRESENT(pw_env_external)) pw_env => pw_env_external
     522      127039 :       IF (PRESENT(task_list_external)) task_list => task_list_external
     523      127039 :       IF (PRESENT(para_env_external)) para_env => para_env_external
     524             : 
     525      127039 :       nspins = dft_control%nspins
     526      127039 :       nimg = dft_control%nimages
     527      127039 :       gapw = dft_control%qs_control%gapw
     528      127039 :       gapw_xc = dft_control%qs_control%gapw_xc
     529             : 
     530      127039 :       CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     531      279140 :       DO ispin = 1, nspins
     532      152101 :          rho_ao => rho_ao_kp(ispin, :)
     533             :          CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     534             :                                  rho=rho_r(ispin), &
     535             :                                  rho_gspace=rho_g(ispin), &
     536             :                                  total_rho=tot_rho_r(ispin), &
     537             :                                  ks_env=ks_env, soft_valid=gapw, &
     538             :                                  task_list_external=task_list_external, &
     539      279140 :                                  pw_env_external=pw_env_external)
     540             :       END DO
     541      127039 :       CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     542             : 
     543      127039 :       IF (gapw_xc) THEN
     544        3644 :          IF (PRESENT(rho_xc_external)) THEN
     545         606 :             rho_xc => rho_xc_external
     546             :          ELSE
     547        3038 :             CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
     548             :          END IF
     549             :          CALL qs_rho_get(rho_xc, &
     550             :                          rho_ao_kp=rho_xc_ao, &
     551             :                          rho_r=rho_xc_r, &
     552             :                          rho_g=rho_xc_g, &
     553        3644 :                          tot_rho_r=tot_rho_r_xc)
     554             :          ! copy rho_ao into rho_xc_ao
     555        7696 :          DO ispin = 1, nspins
     556       13776 :             DO img = 1, nimg
     557       10132 :                CALL dbcsr_copy(rho_xc_ao(ispin, img)%matrix, rho_ao_kp(ispin, img)%matrix)
     558             :             END DO
     559             :          END DO
     560        7696 :          DO ispin = 1, nspins
     561        4052 :             rho_ao => rho_xc_ao(ispin, :)
     562             :             CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     563             :                                     rho=rho_xc_r(ispin), &
     564             :                                     rho_gspace=rho_xc_g(ispin), &
     565             :                                     total_rho=tot_rho_r_xc(ispin), &
     566             :                                     ks_env=ks_env, soft_valid=gapw_xc, &
     567             :                                     task_list_external=task_list_external_soft, &
     568        7696 :                                     pw_env_external=pw_env_external)
     569             :          END DO
     570        3644 :          CALL qs_rho_set(rho_xc, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     571             :       END IF
     572             : 
     573             :       ! GAPW o GAPW_XC require the calculation of hard and soft local densities
     574      127039 :       IF (gapw .OR. gapw_xc) THEN
     575             :          CALL get_qs_env(qs_env=qs_env, &
     576             :                          rho_atom_set=rho_atom_set, &
     577             :                          qs_kind_set=qs_kind_set, &
     578       20208 :                          oce=oce, sab_orb=sab)
     579       20208 :          IF (PRESENT(local_rho_set)) rho_atom_set => local_rho_set%rho_atom_set
     580       20208 :          CPASSERT(ASSOCIATED(rho_atom_set))
     581       20208 :          CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     582       20208 :          CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp, rho_atom_set, qs_kind_set, oce, sab, para_env)
     583             :       END IF
     584             : 
     585      127039 :       IF (.NOT. gapw_xc) THEN
     586             :          ! if needed compute also the gradient of the density
     587      123395 :          IF (dft_control%drho_by_collocation) THEN
     588           0 :             CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     589           0 :             CPASSERT(.NOT. PRESENT(task_list_external))
     590           0 :             DO ispin = 1, nspins
     591           0 :                rho_ao => rho_ao_kp(ispin, :)
     592             :                CALL calculate_drho_elec(matrix_p_kp=rho_ao, &
     593             :                                         drho=drho_r(:, ispin), &
     594             :                                         drho_gspace=drho_g(:, ispin), &
     595           0 :                                         qs_env=qs_env, soft_valid=gapw)
     596             :             END DO
     597           0 :             CALL qs_rho_set(rho_struct, drho_r_valid=.TRUE., drho_g_valid=.TRUE.)
     598             :          END IF
     599             :          ! if needed compute also the kinetic energy density
     600      123395 :          IF (dft_control%use_kinetic_energy_density) THEN
     601        3428 :             CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     602        7374 :             DO ispin = 1, nspins
     603        3946 :                rho_ao => rho_ao_kp(ispin, :)
     604             :                CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     605             :                                        rho=tau_r(ispin), &
     606             :                                        rho_gspace=tau_g(ispin), &
     607             :                                        total_rho=dum, & ! presumably not meaningful
     608             :                                        ks_env=ks_env, soft_valid=gapw, &
     609             :                                        compute_tau=.TRUE., &
     610             :                                        task_list_external=task_list_external, &
     611        7374 :                                        pw_env_external=pw_env_external)
     612             :             END DO
     613        3428 :             CALL qs_rho_set(rho_struct, tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
     614             :          END IF
     615             :       ELSE
     616             :          CALL qs_rho_get(rho_xc, &
     617             :                          drho_r=drho_xc_r, &
     618             :                          drho_g=drho_xc_g, &
     619             :                          tau_r=tau_xc_r, &
     620        3644 :                          tau_g=tau_xc_g)
     621             :          ! if needed compute also the gradient of the density
     622        3644 :          IF (dft_control%drho_by_collocation) THEN
     623           0 :             CPASSERT(.NOT. PRESENT(task_list_external))
     624           0 :             DO ispin = 1, nspins
     625           0 :                rho_ao => rho_xc_ao(ispin, :)
     626             :                CALL calculate_drho_elec(matrix_p_kp=rho_ao, &
     627             :                                         drho=drho_xc_r(:, ispin), &
     628             :                                         drho_gspace=drho_xc_g(:, ispin), &
     629           0 :                                         qs_env=qs_env, soft_valid=gapw_xc)
     630             :             END DO
     631           0 :             CALL qs_rho_set(rho_xc, drho_r_valid=.TRUE., drho_g_valid=.TRUE.)
     632             :          END IF
     633             :          ! if needed compute also the kinetic energy density
     634        3644 :          IF (dft_control%use_kinetic_energy_density) THEN
     635          72 :             DO ispin = 1, nspins
     636          36 :                rho_ao => rho_xc_ao(ispin, :)
     637             :                CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     638             :                                        rho=tau_xc_r(ispin), &
     639             :                                        rho_gspace=tau_xc_g(ispin), &
     640             :                                        ks_env=ks_env, soft_valid=gapw_xc, &
     641             :                                        compute_tau=.TRUE., &
     642             :                                        task_list_external=task_list_external_soft, &
     643          72 :                                        pw_env_external=pw_env_external)
     644             :             END DO
     645          36 :             CALL qs_rho_set(rho_xc, tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
     646             :          END IF
     647             :       END IF
     648             : 
     649      127039 :       CALL timestop(handle)
     650             : 
     651      127039 :    END SUBROUTINE qs_rho_update_rho_low
     652             : 
     653             : ! **************************************************************************************************
     654             : !> \brief updates rho_r and rho_g to the rho%rho_ao.
     655             : !>      if use_kinetic_energy_density also computes tau_r and tau_g
     656             : !> \param rho_struct the rho structure that should be updated
     657             : !> \param qs_env the qs_env rho_struct refers to
     658             : !>        the integrated charge in r space
     659             : !> \param pw_env_external    external plane wave environment
     660             : !> \param task_list_external external task list
     661             : !> \param para_env_external ...
     662             : !> \param tddfpt_lri_env ...
     663             : !> \param tddfpt_lri_density ...
     664             : !> \par History
     665             : !>      08.2002 created [fawzi]
     666             : !> \author Fawzi Mohamed
     667             : ! **************************************************************************************************
     668         206 :    SUBROUTINE qs_rho_update_tddfpt(rho_struct, qs_env, pw_env_external, task_list_external, &
     669             :                                    para_env_external, tddfpt_lri_env, tddfpt_lri_density)
     670             :       TYPE(qs_rho_type), INTENT(INOUT)                   :: rho_struct
     671             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     672             :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
     673             :       TYPE(task_list_type), OPTIONAL, POINTER            :: task_list_external
     674             :       TYPE(mp_para_env_type), OPTIONAL, POINTER          :: para_env_external
     675             :       TYPE(lri_environment_type), OPTIONAL, POINTER      :: tddfpt_lri_env
     676             :       TYPE(lri_density_type), OPTIONAL, POINTER          :: tddfpt_lri_density
     677             : 
     678             :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_update_tddfpt'
     679             : 
     680             :       INTEGER                                            :: handle, ispin, nspins
     681         206 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     682             :       LOGICAL                                            :: lri_response
     683         206 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_r
     684         206 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     685         206 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
     686         206 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp
     687             :       TYPE(dft_control_type), POINTER                    :: dft_control
     688             :       TYPE(kpoint_type), POINTER                         :: kpoints
     689             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     690         206 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     691             :       TYPE(pw_env_type), POINTER                         :: pw_env
     692         206 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
     693             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     694             :       TYPE(task_list_type), POINTER                      :: task_list
     695             : 
     696         206 :       CALL timeset(routineN, handle)
     697             : 
     698             :       CALL get_qs_env(qs_env, &
     699             :                       ks_env=ks_env, &
     700             :                       dft_control=dft_control, &
     701             :                       atomic_kind_set=atomic_kind_set, &
     702             :                       task_list=task_list, &
     703             :                       para_env=para_env, &
     704         206 :                       pw_env=pw_env)
     705         206 :       IF (PRESENT(pw_env_external)) pw_env => pw_env_external
     706         206 :       IF (PRESENT(task_list_external)) task_list => task_list_external
     707         206 :       IF (PRESENT(para_env_external)) para_env => para_env_external
     708             : 
     709             :       CALL qs_rho_get(rho_struct, &
     710             :                       rho_r=rho_r, &
     711             :                       rho_g=rho_g, &
     712         206 :                       tot_rho_r=tot_rho_r)
     713             : 
     714         206 :       nspins = dft_control%nspins
     715             : 
     716         206 :       lri_response = PRESENT(tddfpt_lri_env)
     717         206 :       IF (lri_response) THEN
     718         206 :          CPASSERT(PRESENT(tddfpt_lri_density))
     719             :       END IF
     720             : 
     721         206 :       CPASSERT(.NOT. dft_control%drho_by_collocation)
     722         206 :       CPASSERT(.NOT. dft_control%use_kinetic_energy_density)
     723         206 :       CPASSERT(.NOT. dft_control%qs_control%gapw)
     724         206 :       CPASSERT(.NOT. dft_control%qs_control%gapw_xc)
     725             : 
     726         206 :       IF (lri_response) THEN
     727         206 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
     728         206 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     729         206 :          CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     730             :          CALL calculate_lri_densities(tddfpt_lri_env, tddfpt_lri_density, qs_env, rho_ao_kp, cell_to_index, &
     731             :                                       lri_rho_struct=rho_struct, &
     732             :                                       atomic_kind_set=atomic_kind_set, &
     733             :                                       para_env=para_env, &
     734         206 :                                       response_density=lri_response)
     735         206 :          CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     736             :       ELSE
     737           0 :          CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     738           0 :          DO ispin = 1, nspins
     739           0 :             rho_ao => rho_ao_kp(ispin, :)
     740             :             CALL calculate_rho_elec(matrix_p_kp=rho_ao, &
     741             :                                     rho=rho_r(ispin), &
     742             :                                     rho_gspace=rho_g(ispin), &
     743             :                                     total_rho=tot_rho_r(ispin), &
     744             :                                     ks_env=ks_env, &
     745             :                                     task_list_external=task_list_external, &
     746           0 :                                     pw_env_external=pw_env_external)
     747             :          END DO
     748           0 :          CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     749             :       END IF
     750             : 
     751         206 :       CALL timestop(handle)
     752             : 
     753         206 :    END SUBROUTINE qs_rho_update_tddfpt
     754             : 
     755             : ! **************************************************************************************************
     756             : !> \brief Allocate a density structure and fill it with data from an input structure
     757             : !>        SIZE(rho_input) == mspin == 1  direct copy
     758             : !>        SIZE(rho_input) == mspin == 2  direct copy of alpha and beta spin
     759             : !>        SIZE(rho_input) == 1 AND mspin == 2  copy rho/2 into alpha and beta spin
     760             : !> \param rho_input ...
     761             : !> \param rho_output ...
     762             : !> \param auxbas_pw_pool ...
     763             : !> \param mspin ...
     764             : ! **************************************************************************************************
     765       25936 :    SUBROUTINE qs_rho_copy(rho_input, rho_output, auxbas_pw_pool, mspin)
     766             : 
     767             :       TYPE(qs_rho_type), INTENT(IN)                      :: rho_input
     768             :       TYPE(qs_rho_type), INTENT(INOUT)                   :: rho_output
     769             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     770             :       INTEGER, INTENT(IN)                                :: mspin
     771             : 
     772             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_rho_copy'
     773             : 
     774             :       INTEGER                                            :: handle, i, j, nspins
     775             :       LOGICAL :: complex_rho_ao, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, rho_r_valid_in, &
     776             :          soft_valid_in, tau_g_valid_in, tau_r_valid_in
     777             :       REAL(KIND=dp)                                      :: ospin
     778       12968 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_g_in, tot_rho_g_out, &
     779       12968 :                                                             tot_rho_r_in, tot_rho_r_out
     780       12968 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
     781       12968 :                                                             rho_ao_out
     782       12968 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp_in
     783       12968 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
     784       12968 :       TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER     :: drho_g_in, drho_g_out
     785       12968 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
     786       12968 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER     :: drho_r_in, drho_r_out
     787             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_r_sccs_in, rho_r_sccs_out
     788             : 
     789       12968 :       CALL timeset(routineN, handle)
     790             : 
     791       12968 :       CPASSERT(mspin == 1 .OR. mspin == 2)
     792       12968 :       ospin = 1._dp/REAL(mspin, KIND=dp)
     793             : 
     794       12968 :       CALL qs_rho_clear(rho_output)
     795             : 
     796       12968 :       NULLIFY (rho_ao_in, rho_ao_kp_in, rho_ao_im_in, rho_r_in, rho_g_in, drho_r_in, &
     797       12968 :                drho_g_in, tau_r_in, tau_g_in, tot_rho_r_in, tot_rho_g_in, rho_r_sccs_in)
     798             : 
     799             :       CALL qs_rho_get(rho_input, &
     800             :                       rho_ao=rho_ao_in, &
     801             :                       rho_ao_kp=rho_ao_kp_in, &
     802             :                       rho_ao_im=rho_ao_im_in, &
     803             :                       rho_r=rho_r_in, &
     804             :                       rho_g=rho_g_in, &
     805             :                       drho_r=drho_r_in, &
     806             :                       drho_g=drho_g_in, &
     807             :                       tau_r=tau_r_in, &
     808             :                       tau_g=tau_g_in, &
     809             :                       tot_rho_r=tot_rho_r_in, &
     810             :                       tot_rho_g=tot_rho_g_in, &
     811             :                       rho_g_valid=rho_g_valid_in, &
     812             :                       rho_r_valid=rho_r_valid_in, &
     813             :                       drho_g_valid=drho_g_valid_in, &
     814             :                       drho_r_valid=drho_r_valid_in, &
     815             :                       tau_r_valid=tau_r_valid_in, &
     816             :                       tau_g_valid=tau_g_valid_in, &
     817             :                       rho_r_sccs=rho_r_sccs_in, &
     818             :                       soft_valid=soft_valid_in, &
     819       12968 :                       complex_rho_ao=complex_rho_ao)
     820             : 
     821       12968 :       NULLIFY (rho_ao_out, rho_ao_im_out, rho_r_out, rho_g_out, drho_r_out, &
     822       12968 :                drho_g_out, tau_r_out, tau_g_out, tot_rho_r_out, tot_rho_g_out, rho_r_sccs_out)
     823             :       ! rho_ao
     824       12968 :       IF (ASSOCIATED(rho_ao_in)) THEN
     825       12968 :          nspins = SIZE(rho_ao_in)
     826       12968 :          CPASSERT(mspin >= nspins)
     827       12968 :          CALL dbcsr_allocate_matrix_set(rho_ao_out, mspin)
     828       12968 :          CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
     829       12968 :          IF (mspin > nspins) THEN
     830        2772 :             DO i = 1, mspin
     831        1848 :                ALLOCATE (rho_ao_out(i)%matrix)
     832        1848 :                CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(1)%matrix, name="RHO copy")
     833        2772 :                CALL dbcsr_scale(rho_ao_out(i)%matrix, ospin)
     834             :             END DO
     835             :          ELSE
     836       26104 :             DO i = 1, nspins
     837       14060 :                ALLOCATE (rho_ao_out(i)%matrix)
     838       26104 :                CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, name="RHO copy")
     839             :             END DO
     840             :          END IF
     841             :       END IF
     842             : 
     843             :       ! rho_ao_kp
     844             :       ! only for non-kp, we could probably just copy this pointer, should work also for non-kp?
     845             :       !IF (ASSOCIATED(rho_ao_kp_in)) THEN
     846             :       !   CPABORT("Copy not available")
     847             :       !END IF
     848             : 
     849             :       ! rho_ao_im
     850       12968 :       IF (ASSOCIATED(rho_ao_im_in)) THEN
     851           0 :          nspins = SIZE(rho_ao_im_in)
     852           0 :          CPASSERT(mspin >= nspins)
     853           0 :          CALL dbcsr_allocate_matrix_set(rho_ao_im_out, mspin)
     854           0 :          CALL qs_rho_set(rho_output, rho_ao_im=rho_ao_im_out)
     855           0 :          IF (mspin > nspins) THEN
     856           0 :             DO i = 1, mspin
     857           0 :                ALLOCATE (rho_ao_im_out(i)%matrix)
     858           0 :                CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(1)%matrix, name="RHO copy")
     859           0 :                CALL dbcsr_scale(rho_ao_im_out(i)%matrix, ospin)
     860             :             END DO
     861             :          ELSE
     862           0 :             DO i = 1, nspins
     863           0 :                ALLOCATE (rho_ao_im_out(i)%matrix)
     864           0 :                CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, name="RHO copy")
     865             :             END DO
     866             :          END IF
     867             :       END IF
     868             : 
     869             :       ! rho_r
     870       12968 :       IF (ASSOCIATED(rho_r_in)) THEN
     871       12968 :          nspins = SIZE(rho_r_in)
     872       12968 :          CPASSERT(mspin >= nspins)
     873       54812 :          ALLOCATE (rho_r_out(mspin))
     874       12968 :          CALL qs_rho_set(rho_output, rho_r=rho_r_out)
     875       12968 :          IF (mspin > nspins) THEN
     876        2772 :             DO i = 1, mspin
     877        1848 :                CALL auxbas_pw_pool%create_pw(rho_r_out(i))
     878        1848 :                CALL pw_copy(rho_r_in(1), rho_r_out(i))
     879        2772 :                CALL pw_scale(rho_r_out(i), ospin)
     880             :             END DO
     881             :          ELSE
     882       26104 :             DO i = 1, nspins
     883       14060 :                CALL auxbas_pw_pool%create_pw(rho_r_out(i))
     884       26104 :                CALL pw_copy(rho_r_in(i), rho_r_out(i))
     885             :             END DO
     886             :          END IF
     887             :       END IF
     888             : 
     889             :       ! rho_g
     890       12968 :       IF (ASSOCIATED(rho_g_in)) THEN
     891       12968 :          nspins = SIZE(rho_g_in)
     892       12968 :          CPASSERT(mspin >= nspins)
     893       54812 :          ALLOCATE (rho_g_out(mspin))
     894       12968 :          CALL qs_rho_set(rho_output, rho_g=rho_g_out)
     895       12968 :          IF (mspin > nspins) THEN
     896        2772 :             DO i = 1, mspin
     897        1848 :                CALL auxbas_pw_pool%create_pw(rho_g_out(i))
     898        1848 :                CALL pw_copy(rho_g_in(1), rho_g_out(i))
     899        2772 :                CALL pw_scale(rho_g_out(i), ospin)
     900             :             END DO
     901             :          ELSE
     902       26104 :             DO i = 1, nspins
     903       14060 :                CALL auxbas_pw_pool%create_pw(rho_g_out(i))
     904       26104 :                CALL pw_copy(rho_g_in(i), rho_g_out(i))
     905             :             END DO
     906             :          END IF
     907             :       END IF
     908             : 
     909             :       ! SCCS
     910       12968 :       IF (ASSOCIATED(rho_r_sccs_in)) THEN
     911           0 :          CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
     912           0 :          CALL auxbas_pw_pool%create_pw(rho_r_sccs_out)
     913           0 :          CALL pw_copy(rho_r_sccs_in, rho_r_sccs_out)
     914             :       END IF
     915             : 
     916             :       ! drho_r
     917       12968 :       IF (ASSOCIATED(drho_r_in)) THEN
     918           0 :          nspins = SIZE(drho_r_in)
     919           0 :          CPASSERT(mspin >= nspins)
     920           0 :          ALLOCATE (drho_r_out(3, mspin))
     921           0 :          CALL qs_rho_set(rho_output, drho_r=drho_r_out)
     922           0 :          IF (mspin > nspins) THEN
     923           0 :             DO j = 1, mspin
     924           0 :                DO i = 1, 3
     925           0 :                   CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
     926           0 :                   CALL pw_copy(drho_r_in(i, 1), drho_r_out(i, j))
     927           0 :                   CALL pw_scale(drho_r_out(i, j), ospin)
     928             :                END DO
     929             :             END DO
     930             :          ELSE
     931           0 :             DO j = 1, nspins
     932           0 :                DO i = 1, 3
     933           0 :                   CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
     934           0 :                   CALL pw_copy(drho_r_in(i, j), drho_r_out(i, j))
     935             :                END DO
     936             :             END DO
     937             :          END IF
     938             :       END IF
     939             : 
     940             :       ! drho_g
     941       12968 :       IF (ASSOCIATED(drho_g_in)) THEN
     942           0 :          nspins = SIZE(drho_g_in)
     943           0 :          CPASSERT(mspin >= nspins)
     944           0 :          ALLOCATE (drho_g_out(3, mspin))
     945           0 :          CALL qs_rho_set(rho_output, drho_g=drho_g_out)
     946           0 :          IF (mspin > nspins) THEN
     947           0 :             DO j = 1, mspin
     948           0 :                DO i = 1, 3
     949           0 :                   CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
     950           0 :                   CALL pw_copy(drho_g_in(i, 1), drho_g_out(i, j))
     951           0 :                   CALL pw_scale(drho_g_out(i, j), ospin)
     952             :                END DO
     953             :             END DO
     954             :          ELSE
     955           0 :             DO j = 1, nspins
     956           0 :                DO i = 1, 3
     957           0 :                   CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
     958           0 :                   CALL pw_copy(drho_g_in(i, j), drho_g_out(i, j))
     959             :                END DO
     960             :             END DO
     961             :          END IF
     962             :       END IF
     963             : 
     964             :       ! tau_r
     965       12968 :       IF (ASSOCIATED(tau_r_in)) THEN
     966           0 :          nspins = SIZE(tau_r_in)
     967           0 :          CPASSERT(mspin >= nspins)
     968           0 :          ALLOCATE (tau_r_out(mspin))
     969           0 :          CALL qs_rho_set(rho_output, tau_r=tau_r_out)
     970           0 :          IF (mspin > nspins) THEN
     971           0 :             DO i = 1, mspin
     972           0 :                CALL auxbas_pw_pool%create_pw(tau_r_out(i))
     973           0 :                CALL pw_copy(tau_r_in(1), tau_r_out(i))
     974           0 :                CALL pw_scale(tau_r_out(i), ospin)
     975             :             END DO
     976             :          ELSE
     977           0 :             DO i = 1, nspins
     978           0 :                CALL auxbas_pw_pool%create_pw(tau_r_out(i))
     979           0 :                CALL pw_copy(tau_r_in(i), tau_r_out(i))
     980             :             END DO
     981             :          END IF
     982             :       END IF
     983             : 
     984             :       ! tau_g
     985       12968 :       IF (ASSOCIATED(tau_g_in)) THEN
     986           0 :          nspins = SIZE(tau_g_in)
     987           0 :          CPASSERT(mspin >= nspins)
     988           0 :          ALLOCATE (tau_g_out(mspin))
     989           0 :          CALL qs_rho_set(rho_output, tau_g=tau_g_out)
     990           0 :          IF (mspin > nspins) THEN
     991           0 :             DO i = 1, mspin
     992           0 :                CALL auxbas_pw_pool%create_pw(tau_g_out(i))
     993           0 :                CALL pw_copy(tau_g_in(1), tau_g_out(i))
     994           0 :                CALL pw_scale(tau_g_out(i), ospin)
     995             :             END DO
     996             :          ELSE
     997           0 :             DO i = 1, nspins
     998           0 :                CALL auxbas_pw_pool%create_pw(tau_g_out(i))
     999           0 :                CALL pw_copy(tau_g_in(i), tau_g_out(i))
    1000             :             END DO
    1001             :          END IF
    1002             :       END IF
    1003             : 
    1004             :       ! tot_rho_r
    1005       12968 :       IF (ASSOCIATED(tot_rho_r_in)) THEN
    1006       12968 :          nspins = SIZE(tot_rho_r_in)
    1007       12968 :          CPASSERT(mspin >= nspins)
    1008       38904 :          ALLOCATE (tot_rho_r_out(mspin))
    1009       12968 :          CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
    1010       12968 :          IF (mspin > nspins) THEN
    1011        2772 :             DO i = 1, mspin
    1012        2772 :                tot_rho_r_out(i) = tot_rho_r_in(1)*ospin
    1013             :             END DO
    1014             :          ELSE
    1015       26104 :             DO i = 1, nspins
    1016       26104 :                tot_rho_r_out(i) = tot_rho_r_in(i)
    1017             :             END DO
    1018             :          END IF
    1019             :       END IF
    1020             : 
    1021             :       ! tot_rho_g
    1022       12968 :       IF (ASSOCIATED(tot_rho_g_in)) THEN
    1023           0 :          nspins = SIZE(tot_rho_g_in)
    1024           0 :          CPASSERT(mspin >= nspins)
    1025           0 :          ALLOCATE (tot_rho_g_out(mspin))
    1026           0 :          CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
    1027           0 :          IF (mspin > nspins) THEN
    1028           0 :             DO i = 1, mspin
    1029           0 :                tot_rho_g_out(i) = tot_rho_g_in(1)*ospin
    1030             :             END DO
    1031             :          ELSE
    1032           0 :             DO i = 1, nspins
    1033           0 :                tot_rho_g_out(i) = tot_rho_g_in(i)
    1034             :             END DO
    1035             :          END IF
    1036             :       END IF
    1037             : 
    1038             :       CALL qs_rho_set(rho_output, &
    1039             :                       rho_g_valid=rho_g_valid_in, &
    1040             :                       rho_r_valid=rho_r_valid_in, &
    1041             :                       drho_g_valid=drho_g_valid_in, &
    1042             :                       drho_r_valid=drho_r_valid_in, &
    1043             :                       tau_r_valid=tau_r_valid_in, &
    1044             :                       tau_g_valid=tau_g_valid_in, &
    1045             :                       soft_valid=soft_valid_in, &
    1046       12968 :                       complex_rho_ao=complex_rho_ao)
    1047             : 
    1048       12968 :       CALL timestop(handle)
    1049             : 
    1050       12968 :    END SUBROUTINE qs_rho_copy
    1051             : 
    1052             : ! **************************************************************************************************
    1053             : !> \brief rhoa = alpha*rhoa+beta*rhob
    1054             : !> \param rhoa ...
    1055             : !> \param rhob ...
    1056             : !> \param alpha ...
    1057             : !> \param beta ...
    1058             : ! **************************************************************************************************
    1059       12936 :    SUBROUTINE qs_rho_scale_and_add(rhoa, rhob, alpha, beta)
    1060             : 
    1061             :       TYPE(qs_rho_type), INTENT(IN)                      :: rhoa, rhob
    1062             :       REAL(KIND=dp), INTENT(IN)                          :: alpha, beta
    1063             : 
    1064             :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_scale_and_add'
    1065             : 
    1066             :       INTEGER                                            :: handle, i, j, nspina, nspinb, nspins
    1067       12936 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_g_a, tot_rho_g_b, tot_rho_r_a, &
    1068       12936 :                                                             tot_rho_r_b
    1069       12936 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao_a, rho_ao_b, rho_ao_im_a, &
    1070       12936 :                                                             rho_ao_im_b
    1071       12936 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g_a, rho_g_b, tau_g_a, tau_g_b
    1072       12936 :       TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER     :: drho_g_a, drho_g_b
    1073       12936 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r_a, rho_r_b, tau_r_a, tau_r_b
    1074       12936 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER     :: drho_r_a, drho_r_b
    1075             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_r_sccs_a, rho_r_sccs_b
    1076             : 
    1077       12936 :       CALL timeset(routineN, handle)
    1078             : 
    1079       12936 :       NULLIFY (rho_ao_a, rho_ao_im_a, rho_r_a, rho_g_a, drho_r_a, &
    1080       12936 :                drho_g_a, tau_r_a, tau_g_a, tot_rho_r_a, tot_rho_g_a, rho_r_sccs_a)
    1081             : 
    1082             :       CALL qs_rho_get(rhoa, &
    1083             :                       rho_ao=rho_ao_a, &
    1084             :                       rho_ao_im=rho_ao_im_a, &
    1085             :                       rho_r=rho_r_a, &
    1086             :                       rho_g=rho_g_a, &
    1087             :                       drho_r=drho_r_a, &
    1088             :                       drho_g=drho_g_a, &
    1089             :                       tau_r=tau_r_a, &
    1090             :                       tau_g=tau_g_a, &
    1091             :                       tot_rho_r=tot_rho_r_a, &
    1092             :                       tot_rho_g=tot_rho_g_a, &
    1093       12936 :                       rho_r_sccs=rho_r_sccs_a)
    1094             : 
    1095       12936 :       NULLIFY (rho_ao_b, rho_ao_im_b, rho_r_b, rho_g_b, drho_r_b, &
    1096       12936 :                drho_g_b, tau_r_b, tau_g_b, tot_rho_r_b, tot_rho_g_b, rho_r_sccs_b)
    1097             : 
    1098             :       CALL qs_rho_get(rhob, &
    1099             :                       rho_ao=rho_ao_b, &
    1100             :                       rho_ao_im=rho_ao_im_b, &
    1101             :                       rho_r=rho_r_b, &
    1102             :                       rho_g=rho_g_b, &
    1103             :                       drho_r=drho_r_b, &
    1104             :                       drho_g=drho_g_b, &
    1105             :                       tau_r=tau_r_b, &
    1106             :                       tau_g=tau_g_b, &
    1107             :                       tot_rho_r=tot_rho_r_b, &
    1108             :                       tot_rho_g=tot_rho_g_b, &
    1109       12936 :                       rho_r_sccs=rho_r_sccs_b)
    1110             :       ! rho_ao
    1111       12936 :       IF (ASSOCIATED(rho_ao_a) .AND. ASSOCIATED(rho_ao_b)) THEN
    1112       12936 :          nspina = SIZE(rho_ao_a)
    1113       12936 :          nspinb = SIZE(rho_ao_b)
    1114       12936 :          nspins = MIN(nspina, nspinb)
    1115       27888 :          DO i = 1, nspins
    1116       27888 :             CALL dbcsr_add(rho_ao_a(i)%matrix, rho_ao_b(i)%matrix, alpha, beta)
    1117             :          END DO
    1118             :       END IF
    1119             : 
    1120             :       ! rho_ao_im
    1121       12936 :       IF (ASSOCIATED(rho_ao_im_a) .AND. ASSOCIATED(rho_ao_im_b)) THEN
    1122           0 :          nspina = SIZE(rho_ao_im_a)
    1123           0 :          nspinb = SIZE(rho_ao_im_b)
    1124           0 :          nspins = MIN(nspina, nspinb)
    1125           0 :          DO i = 1, nspins
    1126           0 :             CALL dbcsr_add(rho_ao_im_a(i)%matrix, rho_ao_im_b(i)%matrix, alpha, beta)
    1127             :          END DO
    1128             :       END IF
    1129             : 
    1130             :       ! rho_r
    1131       12936 :       IF (ASSOCIATED(rho_r_a) .AND. ASSOCIATED(rho_r_b)) THEN
    1132       12936 :          nspina = SIZE(rho_ao_a)
    1133       12936 :          nspinb = SIZE(rho_ao_b)
    1134       12936 :          nspins = MIN(nspina, nspinb)
    1135       27888 :          DO i = 1, nspins
    1136       27888 :             CALL pw_axpy(rho_r_b(i), rho_r_a(i), beta, alpha)
    1137             :          END DO
    1138             :       END IF
    1139             : 
    1140             :       ! rho_g
    1141       12936 :       IF (ASSOCIATED(rho_g_a) .AND. ASSOCIATED(rho_g_b)) THEN
    1142       12936 :          nspina = SIZE(rho_ao_a)
    1143       12936 :          nspinb = SIZE(rho_ao_b)
    1144       12936 :          nspins = MIN(nspina, nspinb)
    1145       27888 :          DO i = 1, nspins
    1146       27888 :             CALL pw_axpy(rho_g_b(i), rho_g_a(i), beta, alpha)
    1147             :          END DO
    1148             :       END IF
    1149             : 
    1150             :       ! SCCS
    1151       12936 :       IF (ASSOCIATED(rho_r_sccs_a) .AND. ASSOCIATED(rho_r_sccs_b)) THEN
    1152           0 :          CALL pw_axpy(rho_r_sccs_b, rho_r_sccs_a, beta, alpha)
    1153             :       END IF
    1154             : 
    1155             :       ! drho_r
    1156       12936 :       IF (ASSOCIATED(drho_r_a) .AND. ASSOCIATED(drho_r_b)) THEN
    1157           0 :          CPASSERT(ALL(SHAPE(drho_r_a) == SHAPE(drho_r_b))) ! not implemented
    1158           0 :          DO j = 1, SIZE(drho_r_a, 2)
    1159           0 :             DO i = 1, SIZE(drho_r_a, 1)
    1160           0 :                CALL pw_axpy(drho_r_b(i, j), drho_r_a(i, j), beta, alpha)
    1161             :             END DO
    1162             :          END DO
    1163             :       END IF
    1164             : 
    1165             :       ! drho_g
    1166       12936 :       IF (ASSOCIATED(drho_g_a) .AND. ASSOCIATED(drho_g_b)) THEN
    1167           0 :          CPASSERT(ALL(SHAPE(drho_g_a) == SHAPE(drho_g_b))) ! not implemented
    1168           0 :          DO j = 1, SIZE(drho_g_a, 2)
    1169           0 :             DO i = 1, SIZE(drho_g_a, 1)
    1170           0 :                CALL pw_axpy(drho_g_b(i, j), drho_g_a(i, j), beta, alpha)
    1171             :             END DO
    1172             :          END DO
    1173             :       END IF
    1174             : 
    1175             :       ! tau_r
    1176       12936 :       IF (ASSOCIATED(tau_r_a) .AND. ASSOCIATED(tau_r_b)) THEN
    1177           0 :          nspina = SIZE(rho_ao_a)
    1178           0 :          nspinb = SIZE(rho_ao_b)
    1179           0 :          nspins = MIN(nspina, nspinb)
    1180           0 :          DO i = 1, nspins
    1181           0 :             CALL pw_axpy(tau_r_b(i), tau_r_a(i), beta, alpha)
    1182             :          END DO
    1183             :       END IF
    1184             : 
    1185             :       ! tau_g
    1186       12936 :       IF (ASSOCIATED(tau_g_a) .AND. ASSOCIATED(tau_g_b)) THEN
    1187           0 :          nspina = SIZE(rho_ao_a)
    1188           0 :          nspinb = SIZE(rho_ao_b)
    1189           0 :          nspins = MIN(nspina, nspinb)
    1190           0 :          DO i = 1, nspins
    1191           0 :             CALL pw_axpy(tau_g_b(i), tau_g_a(i), beta, alpha)
    1192             :          END DO
    1193             :       END IF
    1194             : 
    1195             :       ! tot_rho_r
    1196       12936 :       IF (ASSOCIATED(tot_rho_r_a) .AND. ASSOCIATED(tot_rho_r_b)) THEN
    1197         360 :          nspina = SIZE(rho_ao_a)
    1198         360 :          nspinb = SIZE(rho_ao_b)
    1199         360 :          nspins = MIN(nspina, nspinb)
    1200        1008 :          DO i = 1, nspins
    1201        1008 :             tot_rho_r_a(i) = alpha*tot_rho_r_a(i) + beta*tot_rho_r_b(i)
    1202             :          END DO
    1203             :       END IF
    1204             : 
    1205             :       ! tot_rho_g
    1206       12936 :       IF (ASSOCIATED(tot_rho_g_a) .AND. ASSOCIATED(tot_rho_g_b)) THEN
    1207           0 :          nspina = SIZE(rho_ao_a)
    1208           0 :          nspinb = SIZE(rho_ao_b)
    1209           0 :          nspins = MIN(nspina, nspinb)
    1210           0 :          DO i = 1, nspins
    1211           0 :             tot_rho_g_a(i) = alpha*tot_rho_g_a(i) + beta*tot_rho_g_b(i)
    1212             :          END DO
    1213             :       END IF
    1214             : 
    1215       12936 :       CALL timestop(handle)
    1216             : 
    1217       12936 :    END SUBROUTINE qs_rho_scale_and_add
    1218             : 
    1219             : ! **************************************************************************************************
    1220             : !> \brief Duplicates a pointer physically
    1221             : !> \param rho_input The rho structure to be duplicated
    1222             : !> \param rho_output The duplicate rho structure
    1223             : !> \param qs_env The QS environment from which the auxiliary PW basis-set
    1224             : !>                pool is taken
    1225             : !> \par History
    1226             : !>      07.2005 initial create [tdk]
    1227             : !> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch)
    1228             : !> \note
    1229             : !>      Associated pointers are deallocated, nullified pointers are NOT accepted!
    1230             : ! **************************************************************************************************
    1231           8 :    SUBROUTINE duplicate_rho_type(rho_input, rho_output, qs_env)
    1232             : 
    1233             :       TYPE(qs_rho_type), INTENT(INOUT)                   :: rho_input, rho_output
    1234             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1235             : 
    1236             :       CHARACTER(len=*), PARAMETER :: routineN = 'duplicate_rho_type'
    1237             : 
    1238             :       INTEGER                                            :: handle, i, j, nspins
    1239             :       LOGICAL :: complex_rho_ao_in, drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, &
    1240             :          rho_r_valid_in, soft_valid_in, tau_g_valid_in, tau_r_valid_in
    1241           4 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_g_in, tot_rho_g_out, &
    1242           4 :                                                             tot_rho_r_in, tot_rho_r_out
    1243           4 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao_im_in, rho_ao_im_out, rho_ao_in, &
    1244           4 :                                                             rho_ao_out
    1245             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1246           4 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g_in, rho_g_out, tau_g_in, tau_g_out
    1247           4 :       TYPE(pw_c1d_gs_type), DIMENSION(:, :), POINTER     :: drho_g_in, drho_g_out
    1248             :       TYPE(pw_env_type), POINTER                         :: pw_env
    1249             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1250           4 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r_in, rho_r_out, tau_r_in, tau_r_out
    1251           4 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), POINTER     :: drho_r_in, drho_r_out
    1252             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_r_sccs_in, rho_r_sccs_out
    1253             : 
    1254           4 :       CALL timeset(routineN, handle)
    1255             : 
    1256           4 :       NULLIFY (dft_control, pw_env, auxbas_pw_pool)
    1257           4 :       NULLIFY (rho_ao_in, rho_ao_out, rho_ao_im_in, rho_ao_im_out)
    1258           4 :       NULLIFY (rho_r_in, rho_r_out, rho_g_in, rho_g_out, drho_r_in, drho_r_out)
    1259           4 :       NULLIFY (drho_g_in, drho_g_out, tau_r_in, tau_r_out, tau_g_in, tau_g_out)
    1260           4 :       NULLIFY (tot_rho_r_in, tot_rho_r_out, tot_rho_g_in, tot_rho_g_out)
    1261           4 :       NULLIFY (rho_r_sccs_in, rho_r_sccs_out)
    1262             : 
    1263           4 :       CPASSERT(ASSOCIATED(qs_env))
    1264             : 
    1265           4 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, dft_control=dft_control)
    1266           4 :       CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
    1267           4 :       nspins = dft_control%nspins
    1268             : 
    1269           4 :       CALL qs_rho_clear(rho_output)
    1270             : 
    1271             :       CALL qs_rho_get(rho_input, &
    1272             :                       rho_ao=rho_ao_in, &
    1273             :                       rho_ao_im=rho_ao_im_in, &
    1274             :                       rho_r=rho_r_in, &
    1275             :                       rho_g=rho_g_in, &
    1276             :                       drho_r=drho_r_in, &
    1277             :                       drho_g=drho_g_in, &
    1278             :                       tau_r=tau_r_in, &
    1279             :                       tau_g=tau_g_in, &
    1280             :                       tot_rho_r=tot_rho_r_in, &
    1281             :                       tot_rho_g=tot_rho_g_in, &
    1282             :                       rho_g_valid=rho_g_valid_in, &
    1283             :                       rho_r_valid=rho_r_valid_in, &
    1284             :                       drho_g_valid=drho_g_valid_in, &
    1285             :                       drho_r_valid=drho_r_valid_in, &
    1286             :                       tau_r_valid=tau_r_valid_in, &
    1287             :                       tau_g_valid=tau_g_valid_in, &
    1288             :                       rho_r_sccs=rho_r_sccs_in, &
    1289             :                       soft_valid=soft_valid_in, &
    1290           4 :                       complex_rho_ao=complex_rho_ao_in)
    1291             : 
    1292             :       ! rho_ao
    1293           4 :       IF (ASSOCIATED(rho_ao_in)) THEN
    1294           4 :          CALL dbcsr_allocate_matrix_set(rho_ao_out, nspins)
    1295           4 :          CALL qs_rho_set(rho_output, rho_ao=rho_ao_out)
    1296           8 :          DO i = 1, nspins
    1297           4 :             ALLOCATE (rho_ao_out(i)%matrix)
    1298             :             CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, &
    1299           4 :                             name="myDensityMatrix_for_Spin_"//TRIM(ADJUSTL(cp_to_string(i))))
    1300           8 :             CALL dbcsr_set(rho_ao_out(i)%matrix, 0.0_dp)
    1301             :          END DO
    1302             :       END IF
    1303             : 
    1304             :       ! rho_ao_im
    1305           4 :       IF (ASSOCIATED(rho_ao_im_in)) THEN
    1306           0 :          CALL dbcsr_allocate_matrix_set(rho_ao_im_out, nspins)
    1307           0 :          CALL qs_rho_set(rho_output, rho_ao=rho_ao_im_out)
    1308           0 :          DO i = 1, nspins
    1309           0 :             ALLOCATE (rho_ao_im_out(i)%matrix)
    1310             :             CALL dbcsr_copy(rho_ao_im_out(i)%matrix, rho_ao_im_in(i)%matrix, &
    1311           0 :                             name="myImagDensityMatrix_for_Spin_"//TRIM(ADJUSTL(cp_to_string(i))))
    1312           0 :             CALL dbcsr_set(rho_ao_im_out(i)%matrix, 0.0_dp)
    1313             :          END DO
    1314             :       END IF
    1315             : 
    1316             :       ! rho_r
    1317           4 :       IF (ASSOCIATED(rho_r_in)) THEN
    1318          16 :          ALLOCATE (rho_r_out(nspins))
    1319           4 :          CALL qs_rho_set(rho_output, rho_r=rho_r_out)
    1320           8 :          DO i = 1, nspins
    1321           4 :             CALL auxbas_pw_pool%create_pw(rho_r_out(i))
    1322           8 :             CALL pw_copy(rho_r_in(i), rho_r_out(i))
    1323             :          END DO
    1324             :       END IF
    1325             : 
    1326             :       ! rho_g
    1327           4 :       IF (ASSOCIATED(rho_g_in)) THEN
    1328          16 :          ALLOCATE (rho_g_out(nspins))
    1329           4 :          CALL qs_rho_set(rho_output, rho_g=rho_g_out)
    1330           8 :          DO i = 1, nspins
    1331           4 :             CALL auxbas_pw_pool%create_pw(rho_g_out(i))
    1332           8 :             CALL pw_copy(rho_g_in(i), rho_g_out(i))
    1333             :          END DO
    1334             :       END IF
    1335             : 
    1336             :       ! SCCS
    1337           4 :       IF (ASSOCIATED(rho_r_sccs_in)) THEN
    1338           0 :          CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out)
    1339           0 :          CALL auxbas_pw_pool%create_pw(rho_r_sccs_out)
    1340           0 :          CALL pw_copy(rho_r_sccs_in, rho_r_sccs_out)
    1341             :       END IF
    1342             : 
    1343             :       ! drho_r and drho_g are only needed if calculated by collocation
    1344           4 :       IF (dft_control%drho_by_collocation) THEN
    1345             :          ! drho_r
    1346           0 :          IF (ASSOCIATED(drho_r_in)) THEN
    1347           0 :             ALLOCATE (drho_r_out(3, nspins))
    1348           0 :             CALL qs_rho_set(rho_output, drho_r=drho_r_out)
    1349           0 :             DO j = 1, nspins
    1350           0 :                DO i = 1, 3
    1351           0 :                   CALL auxbas_pw_pool%create_pw(drho_r_out(i, j))
    1352           0 :                   CALL pw_copy(drho_r_in(i, j), drho_r_out(i, j))
    1353             :                END DO
    1354             :             END DO
    1355             :          END IF
    1356             : 
    1357             :          ! drho_g
    1358           0 :          IF (ASSOCIATED(drho_g_in)) THEN
    1359           0 :             ALLOCATE (drho_g_out(3, nspins))
    1360           0 :             CALL qs_rho_set(rho_output, drho_g=drho_g_out)
    1361           0 :             DO j = 1, nspins
    1362           0 :                DO i = 1, 3
    1363           0 :                   CALL auxbas_pw_pool%create_pw(drho_g_out(i, j))
    1364           0 :                   CALL pw_copy(drho_g_in(i, j), drho_g_out(i, j))
    1365             :                END DO
    1366             :             END DO
    1367             :          END IF
    1368             :       END IF
    1369             : 
    1370             :       ! tau_r and tau_g are only needed in the case of Meta-GGA XC-functionals
    1371             :       ! are used. Therefore they are only allocated if
    1372             :       ! dft_control%use_kinetic_energy_density is true
    1373           4 :       IF (dft_control%use_kinetic_energy_density) THEN
    1374             :          ! tau_r
    1375           0 :          IF (ASSOCIATED(tau_r_in)) THEN
    1376           0 :             ALLOCATE (tau_r_out(nspins))
    1377           0 :             CALL qs_rho_set(rho_output, tau_r=tau_r_out)
    1378           0 :             DO i = 1, nspins
    1379           0 :                CALL auxbas_pw_pool%create_pw(tau_r_out(i))
    1380           0 :                CALL pw_copy(tau_r_in(i), tau_r_out(i))
    1381             :             END DO
    1382             :          END IF
    1383             : 
    1384             :          ! tau_g
    1385           0 :          IF (ASSOCIATED(tau_g_in)) THEN
    1386           0 :             ALLOCATE (tau_g_out(nspins))
    1387           0 :             CALL qs_rho_set(rho_output, tau_g=tau_g_out)
    1388           0 :             DO i = 1, nspins
    1389           0 :                CALL auxbas_pw_pool%create_pw(tau_g_out(i))
    1390           0 :                CALL pw_copy(tau_g_in(i), tau_g_out(i))
    1391             :             END DO
    1392             :          END IF
    1393             :       END IF
    1394             : 
    1395             :       CALL qs_rho_set(rho_output, &
    1396             :                       rho_g_valid=rho_g_valid_in, &
    1397             :                       rho_r_valid=rho_r_valid_in, &
    1398             :                       drho_g_valid=drho_g_valid_in, &
    1399             :                       drho_r_valid=drho_r_valid_in, &
    1400             :                       tau_r_valid=tau_r_valid_in, &
    1401             :                       tau_g_valid=tau_g_valid_in, &
    1402             :                       soft_valid=soft_valid_in, &
    1403           4 :                       complex_rho_ao=complex_rho_ao_in)
    1404             : 
    1405             :       ! tot_rho_r
    1406           4 :       IF (ASSOCIATED(tot_rho_r_in)) THEN
    1407          12 :          ALLOCATE (tot_rho_r_out(nspins))
    1408           4 :          CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out)
    1409           8 :          DO i = 1, nspins
    1410           8 :             tot_rho_r_out(i) = tot_rho_r_in(i)
    1411             :          END DO
    1412             :       END IF
    1413             : 
    1414             :       ! tot_rho_g
    1415           4 :       IF (ASSOCIATED(tot_rho_g_in)) THEN
    1416           0 :          ALLOCATE (tot_rho_g_out(nspins))
    1417           0 :          CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out)
    1418           0 :          DO i = 1, nspins
    1419           0 :             tot_rho_g_out(i) = tot_rho_g_in(i)
    1420             :          END DO
    1421             : 
    1422             :       END IF
    1423             : 
    1424           4 :       CALL timestop(handle)
    1425             : 
    1426           4 :    END SUBROUTINE duplicate_rho_type
    1427             : 
    1428             : ! **************************************************************************************************
    1429             : !> \brief (Re-)allocates rho_ao_im from real part rho_ao
    1430             : !> \param rho ...
    1431             : !> \param qs_env ...
    1432             : ! **************************************************************************************************
    1433          94 :    SUBROUTINE allocate_rho_ao_imag_from_real(rho, qs_env)
    1434             :       TYPE(qs_rho_type), POINTER                         :: rho
    1435             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1436             : 
    1437             :       CHARACTER(LEN=default_string_length)               :: headline
    1438             :       INTEGER                                            :: i, ic, nimages, nspins
    1439          94 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_im_kp, rho_ao_kp
    1440             :       TYPE(dbcsr_type), POINTER                          :: template
    1441             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1442             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1443          94 :          POINTER                                         :: sab_orb
    1444             : 
    1445          94 :       NULLIFY (rho_ao_im_kp, rho_ao_kp, dft_control, template, sab_orb)
    1446             : 
    1447             :       CALL get_qs_env(qs_env, &
    1448             :                       dft_control=dft_control, &
    1449          94 :                       sab_orb=sab_orb)
    1450             : 
    1451          94 :       CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_im_kp, rho_ao_kp=rho_ao_kp)
    1452             : 
    1453          94 :       nspins = dft_control%nspins
    1454          94 :       nimages = dft_control%nimages
    1455             : 
    1456          94 :       CPASSERT(nspins .EQ. SIZE(rho_ao_kp, 1))
    1457          94 :       CPASSERT(nimages .EQ. SIZE(rho_ao_kp, 2))
    1458             : 
    1459          94 :       CALL dbcsr_allocate_matrix_set(rho_ao_im_kp, nspins, nimages)
    1460          94 :       CALL qs_rho_set(rho, rho_ao_im_kp=rho_ao_im_kp)
    1461         202 :       DO i = 1, nspins
    1462         310 :          DO ic = 1, nimages
    1463         108 :             IF (nspins > 1) THEN
    1464          28 :                IF (i == 1) THEN
    1465          14 :                   headline = "IMAGINARY PART OF DENSITY MATRIX FOR ALPHA SPIN"
    1466             :                ELSE
    1467          14 :                   headline = "IMAGINARY PART OF DENSITY MATRIX FOR BETA SPIN"
    1468             :                END IF
    1469             :             ELSE
    1470          80 :                headline = "IMAGINARY PART OF DENSITY MATRIX"
    1471             :             END IF
    1472         108 :             ALLOCATE (rho_ao_im_kp(i, ic)%matrix)
    1473         108 :             template => rho_ao_kp(i, ic)%matrix ! base on real part, but anti-symmetric
    1474             :             CALL dbcsr_create(matrix=rho_ao_im_kp(i, ic)%matrix, template=template, &
    1475         108 :                               name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
    1476         108 :             CALL cp_dbcsr_alloc_block_from_nbl(rho_ao_im_kp(i, ic)%matrix, sab_orb)
    1477         216 :             CALL dbcsr_set(rho_ao_im_kp(i, ic)%matrix, 0.0_dp)
    1478             :          END DO
    1479             :       END DO
    1480             : 
    1481          94 :    END SUBROUTINE allocate_rho_ao_imag_from_real
    1482             : 
    1483             : END MODULE qs_rho_methods

Generated by: LCOV version 1.15