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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : ! **************************************************************************************************
       8              : MODULE qs_rho0_ggrid
       9              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      10              :                                               get_atomic_kind
      11              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      12              :                                               gto_basis_set_type
      13              :    USE cell_types,                      ONLY: cell_type,&
      14              :                                               pbc
      15              :    USE cp_control_types,                ONLY: dft_control_type
      16              :    USE gaussian_gridlevels,             ONLY: gaussian_gridlevel
      17              :    USE grid_api,                        ONLY: GRID_FUNC_AB,&
      18              :                                               collocate_pgf_product
      19              :    USE kinds,                           ONLY: dp
      20              :    USE memory_utilities,                ONLY: reallocate
      21              :    USE message_passing,                 ONLY: mp_para_env_type
      22              :    USE orbital_pointers,                ONLY: indco,&
      23              :                                               nco,&
      24              :                                               ncoset,&
      25              :                                               nso,&
      26              :                                               nsoset
      27              :    USE orbital_transformation_matrices, ONLY: orbtramat
      28              :    USE particle_types,                  ONLY: particle_type
      29              :    USE pw_env_types,                    ONLY: pw_env_get,&
      30              :                                               pw_env_type
      31              :    USE pw_methods,                      ONLY: pw_axpy,&
      32              :                                               pw_copy,&
      33              :                                               pw_integrate_function,&
      34              :                                               pw_scale,&
      35              :                                               pw_transfer,&
      36              :                                               pw_zero
      37              :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      38              :                                               pw_pool_type
      39              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      40              :                                               pw_r3d_rs_type
      41              :    USE qs_environment_types,            ONLY: get_qs_env,&
      42              :                                               qs_environment_type
      43              :    USE qs_force_types,                  ONLY: qs_force_type
      44              :    USE qs_harmonics_atom,               ONLY: get_none0_cg_list,&
      45              :                                               harmonics_atom_type
      46              :    USE qs_integrate_potential,          ONLY: integrate_pgf_product
      47              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      48              :                                               qs_kind_type
      49              :    USE qs_local_rho_types,              ONLY: get_local_rho,&
      50              :                                               local_rho_type
      51              :    USE qs_rho0_types,                   ONLY: get_rho0_mpole,&
      52              :                                               rho0_mpole_type
      53              :    USE qs_rho_atom_types,               ONLY: get_rho_atom,&
      54              :                                               rho_atom_coeff,&
      55              :                                               rho_atom_type
      56              :    USE realspace_grid_types,            ONLY: realspace_grid_desc_p_type,&
      57              :                                               realspace_grid_desc_type,&
      58              :                                               realspace_grid_type,&
      59              :                                               rs_grid_create,&
      60              :                                               rs_grid_release,&
      61              :                                               rs_grid_zero,&
      62              :                                               transfer_pw2rs,&
      63              :                                               transfer_rs2pw
      64              :    USE util,                            ONLY: get_limit
      65              :    USE virial_types,                    ONLY: virial_type
      66              : #include "./base/base_uses.f90"
      67              : 
      68              :    IMPLICIT NONE
      69              : 
      70              :    PRIVATE
      71              : 
      72              :    ! Global parameters (only in this module)
      73              : 
      74              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_ggrid'
      75              : 
      76              :    ! Public subroutines
      77              : 
      78              :    PUBLIC :: put_rho0_on_grid, rho0_s_grid_create, integrate_vhg0_rspace
      79              : 
      80              : CONTAINS
      81              : 
      82              : ! **************************************************************************************************
      83              : !> \brief ...
      84              : !> \param qs_env ...
      85              : !> \param rho0 ...
      86              : !> \param tot_rs_int ...
      87              : !> \param my_pools ...
      88              : !> \param my_rs_grids ...
      89              : !> \param my_rs_descs ...
      90              : ! **************************************************************************************************
      91        16842 :    SUBROUTINE put_rho0_on_grid(qs_env, rho0, tot_rs_int, my_pools, my_rs_grids, my_rs_descs)
      92              : 
      93              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      94              :       TYPE(rho0_mpole_type), POINTER                     :: rho0
      95              :       REAL(KIND=dp), INTENT(OUT)                         :: tot_rs_int
      96              :       TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
      97              :          POINTER                                         :: my_pools
      98              :       TYPE(realspace_grid_type), DIMENSION(:), &
      99              :          OPTIONAL, POINTER                               :: my_rs_grids
     100              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     101              :          OPTIONAL, POINTER                               :: my_rs_descs
     102              : 
     103              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'put_rho0_on_grid'
     104              : 
     105              :       INTEGER                                            :: auxbas_grid, handle, iat, iatom, igrid, &
     106              :                                                             ikind, ithread, j, l0_ikind, lmax0, &
     107              :                                                             nat, nch_ik, nch_max, npme
     108        16842 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, cores
     109              :       LOGICAL                                            :: paw_atom
     110              :       REAL(KIND=dp)                                      :: eps_rho_rspace, rpgf0, zet0
     111              :       REAL(KIND=dp), DIMENSION(3)                        :: ra
     112        16842 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Qlm_c
     113        16842 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: pab
     114        16842 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     115              :       TYPE(cell_type), POINTER                           :: cell
     116              :       TYPE(dft_control_type), POINTER                    :: dft_control
     117              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     118        16842 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     119              :       TYPE(pw_c1d_gs_type)                               :: coeff_gspace
     120              :       TYPE(pw_c1d_gs_type), POINTER                      :: rho0_s_gs
     121              :       TYPE(pw_env_type), POINTER                         :: pw_env
     122        16842 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     123              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     124              :       TYPE(pw_r3d_rs_type)                               :: coeff_rspace, rho0_r_tmp
     125              :       TYPE(pw_r3d_rs_type), POINTER                      :: rho0_s_rs
     126        16842 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     127              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     128        16842 :          POINTER                                         :: descs
     129              :       TYPE(realspace_grid_desc_type), POINTER            :: desc
     130        16842 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: grids
     131              :       TYPE(realspace_grid_type), POINTER                 :: rs_grid
     132              : 
     133        16842 :       CALL timeset(routineN, handle)
     134              : 
     135        16842 :       NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, Qlm_c)
     136              : 
     137        16842 :       NULLIFY (dft_control, pw_env, particle_set, para_env, cell, rho0_s_gs, rho0_s_rs)
     138              :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
     139              :                       particle_set=particle_set, &
     140              :                       atomic_kind_set=atomic_kind_set, &
     141              :                       qs_kind_set=qs_kind_set, &
     142              :                       para_env=para_env, &
     143        16842 :                       pw_env=pw_env, cell=cell)
     144        16842 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     145              : 
     146        16842 :       NULLIFY (descs, pw_pools)
     147        16842 :       CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools)
     148        16842 :       auxbas_grid = pw_env%auxbas_grid
     149              : 
     150        16842 :       NULLIFY (rho0_s_gs, rho0_s_rs)
     151              :       CALL get_rho0_mpole(rho0_mpole=rho0, lmax_0=lmax0, &
     152              :                           zet0_h=zet0, igrid_zet0_s=igrid, &
     153              :                           rho0_s_gs=rho0_s_gs, &
     154        16842 :                           rho0_s_rs=rho0_s_rs)
     155              : 
     156              :       ! *** set up the rs grid at level igrid
     157        16842 :       NULLIFY (rs_grid, desc, pw_pool)
     158              :       ! IF present, overwrite qs grid for new pool
     159        16842 :       IF (PRESENT(my_pools)) THEN
     160         1098 :          desc => my_rs_descs(igrid)%rs_desc
     161         1098 :          rs_grid => my_rs_grids(igrid)
     162         1098 :          pw_pool => my_pools(igrid)%pool
     163              :       ELSE
     164        15744 :          desc => descs(igrid)%rs_desc
     165        15744 :          rs_grid => grids(igrid)
     166        15744 :          pw_pool => pw_pools(igrid)%pool
     167              :       END IF
     168              : 
     169        16842 :       CPASSERT(ASSOCIATED(desc))
     170        16842 :       CPASSERT(ASSOCIATED(pw_pool))
     171              : 
     172        16842 :       IF (igrid /= auxbas_grid) THEN
     173           12 :          CALL pw_pool%create_pw(coeff_rspace)
     174           12 :          CALL pw_pool%create_pw(coeff_gspace)
     175              :       END IF
     176        16842 :       CALL rs_grid_zero(rs_grid)
     177              : 
     178        16842 :       nch_max = ncoset(lmax0)
     179              : 
     180        50526 :       ALLOCATE (pab(nch_max, 1))
     181              : 
     182        50762 :       DO ikind = 1, SIZE(atomic_kind_set)
     183        33920 :          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
     184        33920 :          CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
     185              : 
     186        33920 :          IF (.NOT. paw_atom .AND. dft_control%qs_control%gapw_control%nopaw_as_gpw) CYCLE
     187              : 
     188              :          CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind, &
     189        32286 :                              rpgf0_s=rpgf0)
     190              : 
     191        32286 :          nch_ik = ncoset(l0_ikind)
     192       402100 :          pab = 0.0_dp
     193              : 
     194        32286 :          CALL reallocate(cores, 1, nat)
     195        32286 :          npme = 0
     196        84030 :          cores = 0
     197              : 
     198        84030 :          DO iat = 1, nat
     199        51744 :             iatom = atom_list(iat)
     200        51744 :             ra(:) = pbc(particle_set(iatom)%r, cell)
     201        84030 :             IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN
     202              :                ! replicated realspace grid, split the atoms up between procs
     203        51744 :                IF (MODULO(nat, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN
     204        25872 :                   npme = npme + 1
     205        25872 :                   cores(npme) = iat
     206              :                END IF
     207              :             ELSE
     208            0 :                npme = npme + 1
     209            0 :                cores(npme) = iat
     210              :             END IF
     211              : 
     212              :          END DO
     213              : 
     214              :          ithread = 0
     215       141206 :          DO j = 1, npme
     216              : 
     217        25872 :             iat = cores(j)
     218        25872 :             iatom = atom_list(iat)
     219              : 
     220        25872 :             CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, Qlm_car=Qlm_c)
     221              : 
     222       489292 :             pab(1:nch_ik, 1) = Qlm_c(1:nch_ik)
     223              : 
     224        25872 :             ra(:) = pbc(particle_set(iatom)%r, cell)
     225              : 
     226              :             CALL collocate_pgf_product( &
     227              :                l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
     228              :                ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 1.0_dp, pab, 0, 0, &
     229              :                rs_grid, ga_gb_function=GRID_FUNC_AB, radius=rpgf0, &
     230        59792 :                use_subpatch=.TRUE., subpatch_pattern=0)
     231              : 
     232              :          END DO ! j
     233              : 
     234              :       END DO ! ikind
     235              : 
     236        16842 :       IF (ASSOCIATED(cores)) THEN
     237        16698 :          DEALLOCATE (cores)
     238              :       END IF
     239              : 
     240        16842 :       DEALLOCATE (pab)
     241              : 
     242        16842 :       IF (igrid /= auxbas_grid) THEN
     243           12 :          CALL transfer_rs2pw(rs_grid, coeff_rspace)
     244           12 :          CALL pw_zero(rho0_s_gs)
     245           12 :          CALL pw_transfer(coeff_rspace, coeff_gspace)
     246           12 :          CALL pw_axpy(coeff_gspace, rho0_s_gs)
     247              : 
     248           12 :          tot_rs_int = pw_integrate_function(coeff_rspace, isign=-1)
     249              : 
     250           12 :          CALL pw_pool%give_back_pw(coeff_rspace)
     251           12 :          CALL pw_pool%give_back_pw(coeff_gspace)
     252              :       ELSE
     253              : 
     254        16830 :          CALL pw_pool%create_pw(rho0_r_tmp)
     255              : 
     256        16830 :          CALL transfer_rs2pw(rs_grid, rho0_r_tmp)
     257              : 
     258        16830 :          tot_rs_int = pw_integrate_function(rho0_r_tmp, isign=-1)
     259              : 
     260        16830 :          CALL pw_transfer(rho0_r_tmp, rho0_s_rs)
     261        16830 :          CALL pw_pool%give_back_pw(rho0_r_tmp)
     262              : 
     263        16830 :          CALL pw_zero(rho0_s_gs)
     264        16830 :          CALL pw_transfer(rho0_s_rs, rho0_s_gs)
     265              :       END IF
     266        16842 :       CALL timestop(handle)
     267              : 
     268        33684 :    END SUBROUTINE put_rho0_on_grid
     269              : 
     270              : ! **************************************************************************************************
     271              : !> \brief ...
     272              : !> \param pw_env ...
     273              : !> \param rho0_mpole ...
     274              : ! **************************************************************************************************
     275         2688 :    SUBROUTINE rho0_s_grid_create(pw_env, rho0_mpole)
     276              : 
     277              :       TYPE(pw_env_type), POINTER                         :: pw_env
     278              :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
     279              : 
     280              :       CHARACTER(len=*), PARAMETER :: routineN = 'rho0_s_grid_create'
     281              : 
     282              :       INTEGER                                            :: handle
     283              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     284              : 
     285         2688 :       CALL timeset(routineN, handle)
     286              : 
     287         2688 :       CPASSERT(ASSOCIATED(pw_env))
     288              : 
     289         2688 :       NULLIFY (auxbas_pw_pool)
     290         2688 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     291         2688 :       CPASSERT(ASSOCIATED(auxbas_pw_pool))
     292              : 
     293              :       ! reallocate rho0 on the global grid in real and reciprocal space
     294         2688 :       CPASSERT(ASSOCIATED(rho0_mpole))
     295              : 
     296              :       ! rho0 density in real space
     297         2688 :       IF (ASSOCIATED(rho0_mpole%rho0_s_rs)) THEN
     298          840 :          CALL rho0_mpole%rho0_s_rs%release()
     299              :       ELSE
     300         1848 :          ALLOCATE (rho0_mpole%rho0_s_rs)
     301              :       END IF
     302         2688 :       CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_rs)
     303              : 
     304              :       ! rho0 density in reciprocal space
     305         2688 :       IF (ASSOCIATED(rho0_mpole%rho0_s_gs)) THEN
     306          840 :          CALL rho0_mpole%rho0_s_gs%release()
     307              :       ELSE
     308         1848 :          ALLOCATE (rho0_mpole%rho0_s_gs)
     309              :       END IF
     310         2688 :       CALL auxbas_pw_pool%create_pw(rho0_mpole%rho0_s_gs)
     311              : 
     312              :       ! Find the grid level suitable for rho0_soft
     313         2688 :       rho0_mpole%igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*rho0_mpole%zet0_h)
     314              : 
     315         2688 :       CALL timestop(handle)
     316              : 
     317         2688 :    END SUBROUTINE rho0_s_grid_create
     318              : 
     319              : ! **************************************************************************************************
     320              : !> \brief ...
     321              : !> \param qs_env ...
     322              : !> \param v_rspace ...
     323              : !> \param para_env ...
     324              : !> \param calculate_forces ...
     325              : !> \param local_rho_set ...
     326              : !> \param local_rho_set_2nd ...
     327              : !> \param atener ...
     328              : !> \param kforce ...
     329              : !> \param my_pools ...
     330              : !> \param my_rs_descs ...
     331              : ! **************************************************************************************************
     332        15744 :    SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, &
     333        15744 :                                     local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
     334              : 
     335              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     336              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_rspace
     337              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     338              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     339              :       TYPE(local_rho_type), OPTIONAL, POINTER            :: local_rho_set, local_rho_set_2nd
     340              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atener
     341              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kforce
     342              :       TYPE(pw_pool_p_type), DIMENSION(:), OPTIONAL, &
     343              :          POINTER                                         :: my_pools
     344              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     345              :          OPTIONAL, POINTER                               :: my_rs_descs
     346              : 
     347              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'
     348              : 
     349              :       INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, &
     350              :          ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, lmax0, &
     351              :          lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1, n2, &
     352              :          nat, nch_ik, nch_max, ncurr, nset, nsotot, nspins, num_pe
     353        15744 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: cg_n_list
     354        15744 :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: cg_list
     355        15744 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, lmax, lmin, npgf
     356              :       LOGICAL                                            :: grid_distributed, paw_atom, use_virial
     357              :       REAL(KIND=dp)                                      :: eps_rho_rspace, force_tmp(3), fscale, &
     358              :                                                             ra(3), rpgf0, zet0
     359              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: my_virial_a, my_virial_b
     360        15744 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: hab_sph, norm_l, Qlm
     361        15744 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hab, hdab_sph, intloc, pab
     362        15744 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: a_hdab_sph, hdab, Qlm_gg
     363        15744 :       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: a_hdab
     364        15744 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     365              :       TYPE(cell_type), POINTER                           :: cell
     366              :       TYPE(dft_control_type), POINTER                    :: dft_control
     367              :       TYPE(gto_basis_set_type), POINTER                  :: basis_1c_set
     368              :       TYPE(harmonics_atom_type), POINTER                 :: harmonics
     369        15744 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     370              :       TYPE(pw_c1d_gs_type)                               :: coeff_gaux, coeff_gspace
     371              :       TYPE(pw_env_type), POINTER                         :: pw_env
     372        15744 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     373              :       TYPE(pw_pool_type), POINTER                        :: pw_aux, pw_pool
     374              :       TYPE(pw_r3d_rs_type)                               :: coeff_raux, coeff_rspace
     375        15744 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     376        15744 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     377              :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     378        15744 :          POINTER                                         :: rs_descs
     379              :       TYPE(realspace_grid_desc_type), POINTER            :: rs_desc
     380       314880 :       TYPE(realspace_grid_type)                          :: rs_v
     381              :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
     382        15744 :       TYPE(rho_atom_coeff), DIMENSION(:), POINTER        :: int_local_h, int_local_s
     383        15744 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
     384              :       TYPE(rho_atom_type), POINTER                       :: rho_atom
     385              :       TYPE(virial_type), POINTER                         :: virial
     386              : 
     387        15744 :       CALL timeset(routineN, handle)
     388              : 
     389              :       ! In case of the external density computed forces probably also
     390              :       ! need to be stored outside qs_env. We can then remove the
     391              :       ! attribute 'OPTIONAL' from the argument 'local_rho_set'.
     392              : 
     393              :       ! CPASSERT(.NOT. (calculate_forces .AND. PRESENT(local_rho_set)))
     394              : !      IF (calculate_forces .AND. PRESENT(local_rho_set)) THEN
     395              : !         CPWARN("Forces and External Density!")
     396              : !      END IF
     397              : 
     398        15744 :       NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
     399        15744 :       NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set)
     400              : 
     401              :       CALL get_qs_env(qs_env=qs_env, &
     402              :                       atomic_kind_set=atomic_kind_set, &
     403              :                       qs_kind_set=qs_kind_set, &
     404              :                       cell=cell, &
     405              :                       dft_control=dft_control, &
     406              :                       force=force, pw_env=pw_env, &
     407              :                       rho0_mpole=rho0_mpole, &
     408              :                       rho_atom_set=rho_atom_set, &
     409              :                       particle_set=particle_set, &
     410        15744 :                       virial=virial)
     411              : 
     412        15744 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     413              : 
     414        15744 :       nspins = dft_control%nspins
     415              : 
     416              :       ! The aim of the following code was to return immediately if the subroutine
     417              :       ! was called for triplet excited states in spin-restricted case. This check
     418              :       ! is also performed before invocation of this subroutine. It should be save
     419              :       ! to remove the optional argument 'do_triplet' from the subroutine interface.
     420              :       !my_tddft = PRESENT(local_rho_set)
     421              :       !IF (my_tddft) THEN
     422              :       !   IF (PRESENT(do_triplet)) THEN
     423              :       !      IF (nspins == 1 .AND. do_triplet) RETURN
     424              :       !   ELSE
     425              :       !      IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN
     426              :       !   END IF
     427              :       !END IF
     428              : 
     429        15744 :       IF (PRESENT(local_rho_set)) &
     430         2026 :          CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set)
     431              :       ! Q from rho0_mpole of local_rho_set
     432              :       ! for TDDFT forces we need mixed potential / integral space
     433              :       ! potential stored on local_rho_set_2nd
     434        15744 :       IF (PRESENT(local_rho_set_2nd)) THEN
     435          128 :          CALL get_local_rho(local_rho_set_2nd, rho_atom_set=rho_atom_set)
     436              :       END IF
     437              :       CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, &
     438              :                           zet0_h=zet0, igrid_zet0_s=igrid, &
     439        15744 :                           norm_g0l_h=norm_l)
     440              : 
     441              :       ! Setup of the potential on the multigrids
     442        15744 :       NULLIFY (rs_descs, pw_pools)
     443        15744 :       CPASSERT(ASSOCIATED(pw_env))
     444        15744 :       CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools)
     445              : 
     446              :       ! Assign from pw_env
     447        15744 :       auxbas_grid = pw_env%auxbas_grid
     448              : 
     449              :       ! IF present, overwrite qs grid for new pool
     450              :       ! Get the potential on the right grid
     451        15744 :       IF (PRESENT(my_pools)) THEN
     452           50 :          rs_desc => my_rs_descs(igrid)%rs_desc
     453           50 :          pw_pool => my_pools(igrid)%pool
     454              :       ELSE
     455        15694 :          rs_desc => rs_descs(igrid)%rs_desc
     456        15694 :          pw_pool => pw_pools(igrid)%pool
     457              :       END IF
     458              : 
     459        15744 :       CALL pw_pool%create_pw(coeff_gspace)
     460        15744 :       CALL pw_pool%create_pw(coeff_rspace)
     461              : 
     462        15744 :       IF (igrid /= auxbas_grid) THEN
     463           12 :          pw_aux => pw_pools(auxbas_grid)%pool
     464           12 :          CALL pw_aux%create_pw(coeff_gaux)
     465           12 :          CALL pw_transfer(v_rspace, coeff_gaux)
     466           12 :          CALL pw_copy(coeff_gaux, coeff_gspace)
     467           12 :          CALL pw_transfer(coeff_gspace, coeff_rspace)
     468           12 :          CALL pw_aux%give_back_pw(coeff_gaux)
     469           12 :          CALL pw_aux%create_pw(coeff_raux)
     470           12 :          fscale = coeff_rspace%pw_grid%dvol/coeff_raux%pw_grid%dvol
     471           12 :          CALL pw_scale(coeff_rspace, fscale)
     472           12 :          CALL pw_aux%give_back_pw(coeff_raux)
     473              :       ELSE
     474              : 
     475        15732 :          IF (coeff_gspace%pw_grid%spherical) THEN
     476            0 :             CALL pw_transfer(v_rspace, coeff_gspace)
     477            0 :             CALL pw_transfer(coeff_gspace, coeff_rspace)
     478              :          ELSE
     479        15732 :             CALL pw_copy(v_rspace, coeff_rspace)
     480              :          END IF
     481              :       END IF
     482        15744 :       CALL pw_pool%give_back_pw(coeff_gspace)
     483              : 
     484              :       ! Setup the rs grid at level igrid
     485        15744 :       CALL rs_grid_create(rs_v, rs_desc)
     486        15744 :       CALL rs_grid_zero(rs_v)
     487        15744 :       CALL transfer_pw2rs(rs_v, coeff_rspace)
     488              : 
     489        15744 :       CALL pw_pool%give_back_pw(coeff_rspace)
     490              : 
     491              :       ! Now the potential is on the right grid => integration
     492              : 
     493        15744 :       eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
     494              : 
     495              :       ! Allocate work storage
     496              : 
     497        15744 :       NULLIFY (hab, hab_sph, hdab, hdab_sph, pab, a_hdab, a_hdab_sph)
     498        15744 :       nch_max = ncoset(lmax0)
     499        15744 :       CALL reallocate(hab, 1, nch_max, 1, 1)
     500        15744 :       CALL reallocate(hab_sph, 1, nch_max)
     501        15744 :       CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1)
     502        15744 :       CALL reallocate(hdab_sph, 1, 3, 1, nch_max)
     503        15744 :       CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1)
     504        15744 :       CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max)
     505        15744 :       CALL reallocate(pab, 1, nch_max, 1, 1)
     506              : 
     507        15744 :       ncurr = -1
     508              : 
     509        15744 :       grid_distributed = rs_v%desc%distributed
     510              : 
     511        15744 :       fscale = 1.0_dp
     512        15744 :       IF (PRESENT(kforce)) THEN
     513           40 :          fscale = kforce
     514              :       END IF
     515              : 
     516        47110 :       DO ikind = 1, SIZE(atomic_kind_set, 1)
     517        31366 :          NULLIFY (basis_1c_set, atom_list, harmonics)
     518        31366 :          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
     519              :          CALL get_qs_kind(qs_kind_set(ikind), &
     520              :                           basis_set=basis_1c_set, basis_type="GAPW_1C", &
     521              :                           paw_atom=paw_atom, &
     522        31366 :                           harmonics=harmonics)
     523              : 
     524        31366 :          IF (.NOT. paw_atom) CYCLE
     525              : 
     526        29782 :          NULLIFY (Qlm_gg, lmax, npgf)
     527              :          CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
     528              :                              l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, &  ! Qs different
     529        29782 :                              rpgf0_s=rpgf0)
     530              : 
     531              :          CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
     532              :                                 lmax=lmax, lmin=lmin, &
     533              :                                 maxso=maxso, maxl=maxl, &
     534        29782 :                                 nset=nset, npgf=npgf)
     535              : 
     536        29782 :          nsotot = maxso*nset
     537       119128 :          ALLOCATE (intloc(nsotot, nsotot))
     538              : 
     539              :          ! Initialize the local KS integrals
     540              : 
     541        29782 :          nch_ik = ncoset(l0_ikind)
     542       367852 :          pab = 1.0_dp
     543        29782 :          max_s_harm = harmonics%max_s_harm
     544        29782 :          llmax = harmonics%llmax
     545              : 
     546       178692 :          ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
     547              : 
     548        29782 :          num_pe = para_env%num_pe
     549        29782 :          mepos = para_env%mepos
     550        89346 :          DO j = 0, num_pe - 1
     551        59564 :             bo = get_limit(nat, num_pe, j)
     552        59564 :             IF (.NOT. grid_distributed .AND. j /= mepos) CYCLE
     553              : 
     554        83649 :             DO iat = bo(1), bo(2)
     555        24085 :                iatom = atom_list(iat)
     556        24085 :                ra(:) = pbc(particle_set(iatom)%r, cell)
     557              : 
     558        24085 :                NULLIFY (Qlm)
     559        24085 :                CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, Qlm_tot=Qlm)
     560              : 
     561       293924 :                hab = 0.0_dp
     562      1031186 :                hdab = 0.0_dp
     563     72529203 :                intloc = 0._dp
     564        24085 :                IF (use_virial) THEN
     565          294 :                   my_virial_a = 0.0_dp
     566          294 :                   my_virial_b = 0.0_dp
     567        38808 :                   a_hdab = 0.0_dp
     568              :                END IF
     569              : 
     570              :                CALL integrate_pgf_product( &
     571              :                   l0_ikind, zet0, 0, 0, 0.0_dp, 0, &
     572              :                   ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_v, &
     573              :                   hab, pab, o1=0, o2=0, &
     574              :                   radius=rpgf0, &
     575              :                   calculate_forces=calculate_forces, &
     576              :                   use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, &
     577        24085 :                   hdab=hdab, a_hdab=a_hdab, use_subpatch=.TRUE., subpatch_pattern=0)
     578              : 
     579              :                ! Convert from cartesian to spherical
     580        86892 :                DO lshell = 0, l0_ikind
     581       267657 :                   DO is = 1, nso(lshell)
     582       180765 :                      iso = is + nsoset(lshell - 1)
     583       180765 :                      hab_sph(iso) = 0.0_dp
     584       723060 :                      hdab_sph(1:3, iso) = 0.0_dp
     585      2349945 :                      a_hdab_sph(1:3, 1:3, iso) = 0.0_dp
     586      1059920 :                      DO ic = 1, nco(lshell)
     587       816348 :                         ico = ic + ncoset(lshell - 1)
     588       816348 :                         lx = indco(1, ico)
     589       816348 :                         ly = indco(2, ico)
     590       816348 :                         lz = indco(3, ico)
     591              :                         hab_sph(iso) = hab_sph(iso) + &
     592              :                                        norm_l(lshell)* &
     593              :                                        orbtramat(lshell)%slm(is, ic)* &
     594       816348 :                                        hab(ico, 1)
     595       816348 :                         IF (calculate_forces) THEN
     596              :                            hdab_sph(1:3, iso) = hdab_sph(1:3, iso) + &
     597              :                                                 norm_l(lshell)* &
     598              :                                                 orbtramat(lshell)%slm(is, ic)* &
     599       296112 :                                                 hdab(1:3, ico, 1)
     600              :                         END IF
     601       997113 :                         IF (use_virial) THEN
     602        47040 :                            DO ii = 1, 3
     603       152880 :                            DO i = 1, 3
     604              :                               a_hdab_sph(i, ii, iso) = a_hdab_sph(i, ii, iso) + &
     605              :                                                        norm_l(lshell)* &
     606              :                                                        orbtramat(lshell)%slm(is, ic)* &
     607       141120 :                                                        a_hdab(i, ii, ico, 1)
     608              :                            END DO
     609              :                            END DO
     610              :                         END IF
     611              : 
     612              :                      END DO ! ic
     613              :                   END DO ! is
     614              :                END DO ! lshell
     615              : 
     616              :                m1 = 0
     617        84052 :                DO iset1 = 1, nset
     618              : 
     619              :                   m2 = 0
     620       259120 :                   DO iset2 = 1, nset
     621              :                      CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
     622       199153 :                                             max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
     623       199153 :                      n1 = nsoset(lmax(iset1))
     624       683058 :                      DO ipgf1 = 1, npgf(iset1)
     625       483905 :                         n2 = nsoset(lmax(iset2))
     626      2073385 :                         DO ipgf2 = 1, npgf(iset2)
     627              : 
     628      8786162 :                            DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local)
     629     19704069 :                               DO icg = 1, cg_n_list(iso)
     630     11401812 :                                  iso1 = cg_list(1, icg, iso)
     631     11401812 :                                  iso2 = cg_list(2, icg, iso)
     632              : 
     633     11401812 :                                  ig1 = iso1 + n1*(ipgf1 - 1) + m1
     634     11401812 :                                  ig2 = iso2 + n2*(ipgf2 - 1) + m2
     635              : 
     636     18313742 :                                  intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso) ! potential times Q
     637              : 
     638              :                               END DO ! icg
     639              :                            END DO ! iso
     640              : 
     641              :                         END DO ! ipgf2
     642              :                      END DO ! ipgf1
     643       458273 :                      m2 = m2 + maxso
     644              :                   END DO ! iset2
     645        84052 :                   m1 = m1 + maxso
     646              :                END DO ! iset1
     647              : 
     648        24085 :                IF (grid_distributed) THEN
     649              :                   ! Sum result over all processors
     650            0 :                   CALL para_env%sum(intloc)
     651              :                END IF
     652              : 
     653        24085 :                IF (j == mepos) THEN
     654        24085 :                   rho_atom => rho_atom_set(iatom)
     655        24085 :                   CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_local_h, ga_Vlocal_gb_s=int_local_s)
     656        52426 :                   DO ispin = 1, nspins
     657    163790882 :                      int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
     658    163814967 :                      int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
     659              :                   END DO
     660              :                END IF
     661              : 
     662        24085 :                IF (PRESENT(atener)) THEN
     663          178 :                   DO iso = 1, nsoset(l0_ikind)
     664          178 :                      atener(iatom) = atener(iatom) + 0.5_dp*Qlm(iso)*hab_sph(iso)
     665              :                   END DO
     666              :                END IF
     667              : 
     668        24085 :                IF (calculate_forces) THEN
     669          935 :                   force_tmp(1:3) = 0.0_dp
     670         8854 :                   DO iso = 1, nsoset(l0_ikind)
     671         7919 :                      force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso) ! Q here is from local_rho_set
     672         7919 :                      force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso)
     673         8854 :                      force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso)
     674              :                   END DO
     675         3740 :                   force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + fscale*force_tmp(1:3)
     676              :                END IF
     677        83649 :                IF (use_virial) THEN
     678          294 :                   my_virial_a = 0.0_dp
     679         2940 :                   DO iso = 1, nsoset(l0_ikind)
     680        10878 :                      DO ii = 1, 3
     681        34398 :                      DO i = 1, 3
     682              :                         ! Q from local_rho_set
     683        23814 :                         virial%pv_gapw(i, ii) = virial%pv_gapw(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
     684        31752 :                         virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + fscale*Qlm(iso)*a_hdab_sph(i, ii, iso)
     685              :                      END DO
     686              :                      END DO
     687              :                   END DO
     688              :                END IF
     689              : 
     690              :             END DO
     691              :          END DO
     692              : 
     693        29782 :          DEALLOCATE (intloc)
     694       108258 :          DEALLOCATE (cg_list, cg_n_list)
     695              : 
     696              :       END DO ! ikind
     697              : 
     698        15744 :       CALL rs_grid_release(rs_v)
     699              : 
     700        15744 :       DEALLOCATE (hab, hdab, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph)
     701              : 
     702        15744 :       CALL timestop(handle)
     703              : 
     704        31488 :    END SUBROUTINE integrate_vhg0_rspace
     705              : 
     706              : END MODULE qs_rho0_ggrid
        

Generated by: LCOV version 2.0-1