LCOV - code coverage report
Current view: top level - src - qs_linres_atom_current.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 486 508 95.7 %
Date: 2024-03-28 07:31:50 Functions: 5 5 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 given the response wavefunctions obtained by the application
      10             : !>      of the (rxp), p, and ((dk-dl)xp) operators,
      11             : !>      here the current density vector (jx, jy, jz)
      12             : !>      is computed for the 3 directions of the magnetic field (Bx, By, Bz)
      13             : !> \par History
      14             : !>      created 02-2006 [MI]
      15             : !> \author MI
      16             : ! **************************************************************************************************
      17             : MODULE qs_linres_atom_current
      18             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      19             :                                               get_atomic_kind
      20             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      21             :                                               gto_basis_set_p_type,&
      22             :                                               gto_basis_set_type
      23             :    USE cell_types,                      ONLY: cell_type,&
      24             :                                               pbc
      25             :    USE cp_control_types,                ONLY: dft_control_type
      26             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      27             :    USE dbcsr_api,                       ONLY: dbcsr_get_block_p,&
      28             :                                               dbcsr_p_type
      29             :    USE input_constants,                 ONLY: current_gauge_atom,&
      30             :                                               current_gauge_r,&
      31             :                                               current_gauge_r_and_step_func
      32             :    USE kinds,                           ONLY: dp
      33             :    USE message_passing,                 ONLY: mp_para_env_type
      34             :    USE orbital_pointers,                ONLY: indso,&
      35             :                                               nsoset
      36             :    USE particle_types,                  ONLY: particle_type
      37             :    USE paw_basis_types,                 ONLY: get_paw_basis_info
      38             :    USE qs_environment_types,            ONLY: get_qs_env,&
      39             :                                               qs_environment_type
      40             :    USE qs_grid_atom,                    ONLY: grid_atom_type
      41             :    USE qs_harmonics_atom,               ONLY: get_none0_cg_list,&
      42             :                                               harmonics_atom_type
      43             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      44             :                                               get_qs_kind_set,&
      45             :                                               qs_kind_type
      46             :    USE qs_linres_op,                    ONLY: fac_vecp,&
      47             :                                               set_vecp,&
      48             :                                               set_vecp_rev
      49             :    USE qs_linres_types,                 ONLY: allocate_jrho_atom_rad,&
      50             :                                               allocate_jrho_coeff,&
      51             :                                               current_env_type,&
      52             :                                               get_current_env,&
      53             :                                               jrho_atom_type,&
      54             :                                               set2zero_jrho_atom_rad
      55             :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      56             :                                               neighbor_list_iterate,&
      57             :                                               neighbor_list_iterator_create,&
      58             :                                               neighbor_list_iterator_p_type,&
      59             :                                               neighbor_list_iterator_release,&
      60             :                                               neighbor_list_set_p_type
      61             :    USE qs_oce_methods,                  ONLY: proj_blk
      62             :    USE qs_oce_types,                    ONLY: oce_matrix_type
      63             :    USE qs_rho_atom_types,               ONLY: rho_atom_coeff
      64             :    USE sap_kind_types,                  ONLY: alist_pre_align_blk,&
      65             :                                               alist_type,&
      66             :                                               get_alist
      67             :    USE util,                            ONLY: get_limit
      68             : 
      69             : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
      70             : !$                    omp_get_thread_num, &
      71             : !$                    omp_lock_kind, &
      72             : !$                    omp_init_lock, omp_set_lock, &
      73             : !$                    omp_unset_lock, omp_destroy_lock
      74             : 
      75             : #include "./base/base_uses.f90"
      76             : 
      77             :    IMPLICIT NONE
      78             : 
      79             :    PRIVATE
      80             : 
      81             :    ! *** Public subroutines ***
      82             :    PUBLIC :: calculate_jrho_atom_rad, calculate_jrho_atom, calculate_jrho_atom_coeff
      83             : 
      84             :    ! *** Global parameters (only in this module)
      85             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_atom_current'
      86             : 
      87             : CONTAINS
      88             : 
      89             : ! **************************************************************************************************
      90             : !> \brief Calculate the expansion coefficients for the atomic terms
      91             : !>      of the current densitiy in GAPW
      92             : !> \param qs_env ...
      93             : !> \param current_env ...
      94             : !> \param mat_d0 ...
      95             : !> \param mat_jp ...
      96             : !> \param mat_jp_rii ...
      97             : !> \param mat_jp_riii ...
      98             : !> \param iB ...
      99             : !> \param idir ...
     100             : !> \par History
     101             : !>      07.2006 created [MI]
     102             : !>      02.2009 using new setup of projector-basis overlap [jgh]
     103             : !>      08.2016 add OpenMP [EPCC]
     104             : !>      09.2016 use automatic arrays [M Tucker]
     105             : ! **************************************************************************************************
     106         864 :    SUBROUTINE calculate_jrho_atom_coeff(qs_env, current_env, mat_d0, mat_jp, mat_jp_rii, &
     107             :                                         mat_jp_riii, iB, idir)
     108             : 
     109             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     110             :       TYPE(current_env_type)                             :: current_env
     111             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_d0, mat_jp, mat_jp_rii, mat_jp_riii
     112             :       INTEGER, INTENT(IN)                                :: iB, idir
     113             : 
     114             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_jrho_atom_coeff'
     115             : 
     116             :       INTEGER :: bo(2), handle, iac, iat, iatom, ibc, idir2, ii, iii, ikind, ispin, istat, jatom, &
     117             :          jkind, kac, katom, kbc, kkind, len_CPC, len_PC1, max_gau, max_nsgf, mepos, n_cont_a, &
     118             :          n_cont_b, nat, natom, nkind, nsatbas, nsgfa, nsgfb, nso, nsoctot, nspins, num_pe, &
     119             :          output_unit
     120             :       INTEGER, DIMENSION(3)                              :: cell_b
     121         864 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, list_a, list_b
     122             :       LOGICAL                                            :: den_found, dista, distab, distb, &
     123             :                                                             is_not_associated, paw_atom, &
     124             :                                                             sgf_soft_only_a, sgf_soft_only_b
     125             :       REAL(dp)                                           :: eps_cpc, jmax, nbr_dbl, rab(3), rbc(3)
     126         864 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: a_matrix, b_matrix, c_matrix, d_matrix
     127         864 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: C_coeff_hh_a, C_coeff_hh_b, &
     128         864 :                                                             C_coeff_ss_a, C_coeff_ss_b, r_coef_h, &
     129         864 :                                                             r_coef_s, tmp_coeff, zero_coeff
     130             :       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
     131         864 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     132         864 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_a, mat_b, mat_c, mat_d
     133             :       TYPE(dft_control_type), POINTER                    :: dft_control
     134         864 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     135             :       TYPE(gto_basis_set_type), POINTER                  :: basis_1c_set, basis_set_a, basis_set_b
     136         864 :       TYPE(jrho_atom_type), DIMENSION(:), POINTER        :: jrho1_atom_set
     137             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     138             :       TYPE(neighbor_list_iterator_p_type), &
     139         864 :          DIMENSION(:), POINTER                           :: nl_iterator
     140             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     141         864 :          POINTER                                         :: sab_all
     142             :       TYPE(oce_matrix_type), POINTER                     :: oce
     143         864 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     144         864 :       TYPE(rho_atom_coeff), DIMENSION(:), POINTER        :: a_block, b_block, c_block, d_block, &
     145         864 :                                                             jp2_RARnu, jp_RARnu
     146             : 
     147         864 : !$    INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: proj_blk_lock, alloc_lock
     148             : !$    INTEGER                                            :: lock
     149             : 
     150         864 :       CALL timeset(routineN, handle)
     151             : 
     152         864 :       NULLIFY (atomic_kind_set, qs_kind_set, dft_control, sab_all, jrho1_atom_set, oce, &
     153         864 :                para_env, zero_coeff, tmp_coeff)
     154             : 
     155             :       CALL get_qs_env(qs_env=qs_env, &
     156             :                       atomic_kind_set=atomic_kind_set, &
     157             :                       qs_kind_set=qs_kind_set, &
     158             :                       dft_control=dft_control, &
     159             :                       oce=oce, &
     160             :                       sab_all=sab_all, &
     161         864 :                       para_env=para_env)
     162         864 :       CPASSERT(ASSOCIATED(oce))
     163             : 
     164         864 :       CALL get_current_env(current_env=current_env, jrho1_atom_set=jrho1_atom_set)
     165         864 :       CPASSERT(ASSOCIATED(jrho1_atom_set))
     166             : 
     167             :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
     168         864 :                            maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
     169             : 
     170         864 :       eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
     171             : 
     172         864 :       idir2 = 1
     173         864 :       IF (idir .NE. iB) THEN
     174         576 :          CALL set_vecp_rev(idir, iB, idir2)
     175             :       END IF
     176         864 :       CALL set_vecp(iB, ii, iii)
     177             : 
     178             :       ! Set pointers for the different gauge
     179         864 :       mat_a => mat_d0
     180         864 :       mat_b => mat_jp
     181         864 :       mat_c => mat_jp_rii
     182         864 :       mat_d => mat_jp_riii
     183             : 
     184             :       ! Density-like matrices
     185         864 :       nkind = SIZE(qs_kind_set)
     186         864 :       natom = SIZE(jrho1_atom_set)
     187         864 :       nspins = dft_control%nspins
     188             : 
     189             :       ! Reset CJC coefficients and local density arrays
     190        2466 :       DO ikind = 1, nkind
     191        1602 :          NULLIFY (atom_list)
     192             :          CALL get_atomic_kind(atomic_kind_set(ikind), &
     193             :                               atom_list=atom_list, &
     194        1602 :                               natom=nat)
     195        1602 :          CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
     196             : 
     197             :          ! Quick cycle if needed.
     198        1602 :          IF (.NOT. paw_atom) CYCLE
     199             : 
     200             :          ! Initialize the density matrix-like arrays.
     201        5400 :          DO iat = 1, nat
     202        1692 :             iatom = atom_list(iat)
     203        5940 :             DO ispin = 1, nspins
     204        4338 :                IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)) THEN
     205      726688 :                   jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef = 0.0_dp
     206      726688 :                   jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef = 0.0_dp
     207      726688 :                   jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef = 0.0_dp
     208      726688 :                   jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef = 0.0_dp
     209      726688 :                   jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef = 0.0_dp
     210      726688 :                   jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef = 0.0_dp
     211      726688 :                   jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef = 0.0_dp
     212      726688 :                   jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef = 0.0_dp
     213             :                END IF
     214             :             END DO ! ispin
     215             :          END DO ! iat
     216             :       END DO ! ikind
     217             : 
     218             :       ! Three centers
     219        4194 :       ALLOCATE (basis_set_list(nkind))
     220        2466 :       DO ikind = 1, nkind
     221        1602 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
     222        2466 :          IF (ASSOCIATED(basis_set_a)) THEN
     223        1602 :             basis_set_list(ikind)%gto_basis_set => basis_set_a
     224             :          ELSE
     225           0 :             NULLIFY (basis_set_list(ikind)%gto_basis_set)
     226             :          END IF
     227             :       END DO
     228             : 
     229         864 :       len_PC1 = max_nsgf*max_gau
     230         864 :       len_CPC = max_gau*max_gau
     231             : 
     232             :       num_pe = 1
     233         864 : !$    num_pe = omp_get_max_threads()
     234         864 :       CALL neighbor_list_iterator_create(nl_iterator, sab_all, nthread=num_pe)
     235             : 
     236             : !$OMP PARALLEL DEFAULT( NONE )                                  &
     237             : !$OMP           SHARED( nspins, max_nsgf, max_gau               &
     238             : !$OMP                 , len_PC1, len_CPC                        &
     239             : !$OMP                 , nl_iterator, basis_set_list             &
     240             : !$OMP                 , mat_a, mat_b, mat_c, mat_d              &
     241             : !$OMP                 , nkind, qs_kind_set, eps_cpc, oce        &
     242             : !$OMP                 , ii, iii, jrho1_atom_set                 &
     243             : !$OMP                 , natom, proj_blk_lock, alloc_lock        &
     244             : !$OMP                 )                                         &
     245             : !$OMP          PRIVATE( a_block, b_block, c_block, d_block                      &
     246             : !$OMP                 , jp_RARnu, jp2_RARnu                                     &
     247             : !$OMP                 , a_matrix, b_matrix, c_matrix, d_matrix, istat           &
     248             : !$OMP                 , mepos                                                   &
     249             : !$OMP                 , ikind, jkind, kkind, iatom, jatom, katom                &
     250             : !$OMP                 , cell_b, rab, rbc                                        &
     251             : !$OMP                 , basis_set_a, nsgfa                                      &
     252             : !$OMP                 , basis_set_b, nsgfb                                      &
     253             : !$OMP                 , basis_1c_set, jmax, den_found                          &
     254             : !$OMP                 , nsatbas, nsoctot, nso, paw_atom                         &
     255             : !$OMP                 , iac , alist_ac, kac, n_cont_a, list_a, sgf_soft_only_a  &
     256             : !$OMP                 , ibc , alist_bc, kbc, n_cont_b, list_b, sgf_soft_only_b  &
     257             : !$OMP                 , C_coeff_hh_a, C_coeff_ss_a, dista                       &
     258             : !$OMP                 , C_coeff_hh_b, C_coeff_ss_b, distb                       &
     259             : !$OMP                 , distab                                                  &
     260             : !$OMP                 , r_coef_s, r_coef_h                                      &
     261         864 : !$OMP                 )
     262             : 
     263             :       NULLIFY (a_block, b_block, c_block, d_block)
     264             :       NULLIFY (basis_1c_set, jp_RARnu, jp2_RARnu)
     265             : 
     266             :       ALLOCATE (a_block(nspins), b_block(nspins), c_block(nspins), d_block(nspins), &
     267             :                 jp_RARnu(nspins), jp2_RARnu(nspins), &
     268             :                 STAT=istat)
     269             :       CPASSERT(istat == 0)
     270             : 
     271             :       ALLOCATE (a_matrix(max_nsgf, max_nsgf), b_matrix(max_nsgf, max_nsgf), &
     272             :                 c_matrix(max_nsgf, max_nsgf), d_matrix(max_nsgf, max_nsgf), &
     273             :                 STAT=istat)
     274             :       CPASSERT(istat == 0)
     275             : 
     276             : !$OMP SINGLE
     277             : !$    ALLOCATE (alloc_lock(natom))
     278             : !$    ALLOCATE (proj_blk_lock(nspins*natom))
     279             : !$OMP END SINGLE
     280             : 
     281             : !$OMP DO
     282             : !$    DO lock = 1, natom
     283             : !$       call omp_init_lock(alloc_lock(lock))
     284             : !$    END DO
     285             : !$OMP END DO
     286             : 
     287             : !$OMP DO
     288             : !$    DO lock = 1, nspins*natom
     289             : !$       call omp_init_lock(proj_blk_lock(lock))
     290             : !$    END DO
     291             : !$OMP END DO
     292             : 
     293             :       mepos = 0
     294             : !$    mepos = omp_get_thread_num()
     295             : 
     296             :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     297             : 
     298             :          CALL get_iterator_info(nl_iterator, mepos=mepos, &
     299             :                                 ikind=ikind, jkind=jkind, &
     300             :                                 iatom=iatom, jatom=jatom, cell=cell_b, r=rab)
     301             :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     302             :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     303             :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     304             :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     305             :          nsgfa = basis_set_a%nsgf
     306             :          nsgfb = basis_set_b%nsgf
     307             :          DO ispin = 1, nspins
     308             :             NULLIFY (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef)
     309             :             ALLOCATE (jp_RARnu(ispin)%r_coef(nsgfa, nsgfb), &
     310             :                       jp2_RARnu(ispin)%r_coef(nsgfa, nsgfb))
     311             :          END DO
     312             : 
     313             :          ! Take the block \mu\nu of jpab, jpab_ii and jpab_iii
     314             :          jmax = 0._dp
     315             :          DO ispin = 1, nspins
     316             :             NULLIFY (a_block(ispin)%r_coef)
     317             :             NULLIFY (b_block(ispin)%r_coef)
     318             :             NULLIFY (c_block(ispin)%r_coef)
     319             :             NULLIFY (d_block(ispin)%r_coef)
     320             :             CALL dbcsr_get_block_p(matrix=mat_a(ispin)%matrix, &
     321             :                                    row=iatom, col=jatom, block=a_block(ispin)%r_coef, &
     322             :                                    found=den_found)
     323             :             jmax = jmax + MAXVAL(ABS(a_block(ispin)%r_coef))
     324             :             CALL dbcsr_get_block_p(matrix=mat_b(ispin)%matrix, &
     325             :                                    row=iatom, col=jatom, block=b_block(ispin)%r_coef, &
     326             :                                    found=den_found)
     327             :             jmax = jmax + MAXVAL(ABS(b_block(ispin)%r_coef))
     328             :             CALL dbcsr_get_block_p(matrix=mat_c(ispin)%matrix, &
     329             :                                    row=iatom, col=jatom, block=c_block(ispin)%r_coef, &
     330             :                                    found=den_found)
     331             :             jmax = jmax + MAXVAL(ABS(c_block(ispin)%r_coef))
     332             :             CALL dbcsr_get_block_p(matrix=mat_d(ispin)%matrix, &
     333             :                                    row=iatom, col=jatom, block=d_block(ispin)%r_coef, &
     334             :                                    found=den_found)
     335             :             jmax = jmax + MAXVAL(ABS(d_block(ispin)%r_coef))
     336             :          END DO
     337             : 
     338             :          ! Loop over atoms
     339             :          DO kkind = 1, nkind
     340             :             CALL get_qs_kind(qs_kind_set(kkind), &
     341             :                              basis_set=basis_1c_set, basis_type="GAPW_1C", &
     342             :                              paw_atom=paw_atom)
     343             : 
     344             :             ! Quick cycle if needed.
     345             :             IF (.NOT. paw_atom) CYCLE
     346             : 
     347             :             CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
     348             :             nsoctot = nsatbas
     349             : 
     350             :             iac = ikind + nkind*(kkind - 1)
     351             :             ibc = jkind + nkind*(kkind - 1)
     352             :             IF (.NOT. ASSOCIATED(oce%intac(iac)%alist)) CYCLE
     353             :             IF (.NOT. ASSOCIATED(oce%intac(ibc)%alist)) CYCLE
     354             : 
     355             :             CALL get_alist(oce%intac(iac), alist_ac, iatom)
     356             :             CALL get_alist(oce%intac(ibc), alist_bc, jatom)
     357             :             IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
     358             :             IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
     359             : 
     360             :             DO kac = 1, alist_ac%nclist
     361             :                DO kbc = 1, alist_bc%nclist
     362             :                   IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
     363             : 
     364             :                   IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
     365             :                      IF (jmax*alist_bc%clist(kbc)%maxac*alist_ac%clist(kac)%maxac < eps_cpc) CYCLE
     366             : 
     367             :                      n_cont_a = alist_ac%clist(kac)%nsgf_cnt
     368             :                      n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
     369             :                      sgf_soft_only_a = alist_ac%clist(kac)%sgf_soft_only
     370             :                      sgf_soft_only_b = alist_bc%clist(kbc)%sgf_soft_only
     371             :                      IF (n_cont_a .EQ. 0 .OR. n_cont_b .EQ. 0) CYCLE
     372             : 
     373             :                      ! thanks to the linearity of the response, we
     374             :                      ! can avoid computing soft-soft interactions.
     375             :                      ! those terms are already included in the
     376             :                      ! regular grid.
     377             :                      IF (sgf_soft_only_a .AND. sgf_soft_only_b) CYCLE
     378             : 
     379             :                      list_a => alist_ac%clist(kac)%sgf_list
     380             :                      list_b => alist_bc%clist(kbc)%sgf_list
     381             : 
     382             :                      katom = alist_ac%clist(kac)%catom
     383             : !$                   CALL omp_set_lock(alloc_lock(katom))
     384             :                      IF (.NOT. ASSOCIATED(jrho1_atom_set(katom)%cjc0_h(1)%r_coef)) THEN
     385             :                         CALL allocate_jrho_coeff(jrho1_atom_set, katom, nsoctot)
     386             :                      END IF
     387             : !$                   CALL omp_unset_lock(alloc_lock(katom))
     388             : 
     389             :                      ! Compute the modified Qai matrix as
     390             :                      ! mQai_\mu\nu = Qai_\mu\nu - Qbi_\mu\nu * (R_A-R_\nu)_ii
     391             :                      !             + Qci_\mu\nu * (R_A-R_\nu)_iii
     392             :                      rbc = alist_bc%clist(kbc)%rac
     393             :                      DO ispin = 1, nspins
     394             :                         CALL DCOPY(nsgfa*nsgfb, b_block(ispin)%r_coef(1, 1), 1, &
     395             :                                    jp_RARnu(ispin)%r_coef(1, 1), 1)
     396             :                         CALL DAXPY(nsgfa*nsgfb, -rbc(ii), d_block(ispin)%r_coef(1, 1), 1, &
     397             :                                    jp_RARnu(ispin)%r_coef(1, 1), 1)
     398             :                         CALL DAXPY(nsgfa*nsgfb, rbc(iii), c_block(ispin)%r_coef(1, 1), 1, &
     399             :                                    jp_RARnu(ispin)%r_coef(1, 1), 1)
     400             :                      END DO
     401             : 
     402             :                      ! Get the d_A's for the hard and soft densities.
     403             :                      IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
     404             :                         C_coeff_hh_a => alist_ac%clist(kac)%achint(:, :, 1)
     405             :                         C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
     406             :                         dista = .FALSE.
     407             :                      ELSE
     408             :                         C_coeff_hh_a => alist_ac%clist(kac)%acint(:, :, 1)
     409             :                         C_coeff_ss_a => alist_ac%clist(kac)%acint(:, :, 1)
     410             :                         dista = .TRUE.
     411             :                      END IF
     412             :                      ! Get the d_B's for the hard and soft densities.
     413             :                      IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
     414             :                         C_coeff_hh_b => alist_bc%clist(kbc)%achint(:, :, 1)
     415             :                         C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
     416             :                         distb = .FALSE.
     417             :                      ELSE
     418             :                         C_coeff_hh_b => alist_bc%clist(kbc)%acint(:, :, 1)
     419             :                         C_coeff_ss_b => alist_bc%clist(kbc)%acint(:, :, 1)
     420             :                         distb = .TRUE.
     421             :                      END IF
     422             : 
     423             :                      distab = dista .AND. distb
     424             : 
     425             :                      nso = nsoctot
     426             : 
     427             :                      DO ispin = 1, nspins
     428             : 
     429             :                         ! align the blocks
     430             :                         CALL alist_pre_align_blk(a_block(ispin)%r_coef, SIZE(a_block(ispin)%r_coef, 1), &
     431             :                                                  a_matrix, SIZE(a_matrix, 1), &
     432             :                                                  list_a, n_cont_a, list_b, n_cont_b)
     433             : 
     434             :                         CALL alist_pre_align_blk(jp_RARnu(ispin)%r_coef, SIZE(jp_RARnu(ispin)%r_coef, 1), &
     435             :                                                  b_matrix, SIZE(b_matrix, 1), &
     436             :                                                  list_a, n_cont_a, list_b, n_cont_b)
     437             : 
     438             :                         CALL alist_pre_align_blk(c_block(ispin)%r_coef, SIZE(c_block(ispin)%r_coef, 1), &
     439             :                                                  c_matrix, SIZE(c_matrix, 1), &
     440             :                                                  list_a, n_cont_a, list_b, n_cont_b)
     441             :                         CALL alist_pre_align_blk(d_block(ispin)%r_coef, SIZE(d_block(ispin)%r_coef, 1), &
     442             :                                                  d_matrix, SIZE(d_matrix, 1), &
     443             :                                                  list_a, n_cont_a, list_b, n_cont_b)
     444             : 
     445             : !$                      CALL omp_set_lock(proj_blk_lock((katom - 1)*nspins + ispin))
     446             :                         !------------------------------------------------------------------
     447             :                         ! P_\alpha\alpha'
     448             :                         r_coef_h => jrho1_atom_set(katom)%cjc0_h(ispin)%r_coef
     449             :                         r_coef_s => jrho1_atom_set(katom)%cjc0_s(ispin)%r_coef
     450             :                         CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
     451             :                                       C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
     452             :                                       a_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
     453             :                                       len_PC1, len_CPC, 1.0_dp, distab)
     454             :                         !------------------------------------------------------------------
     455             :                         ! mQai_\alpha\alpha'
     456             :                         r_coef_h => jrho1_atom_set(katom)%cjc_h(ispin)%r_coef
     457             :                         r_coef_s => jrho1_atom_set(katom)%cjc_s(ispin)%r_coef
     458             :                         CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
     459             :                                       C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
     460             :                                       b_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
     461             :                                       len_PC1, len_CPC, 1.0_dp, distab)
     462             :                         !------------------------------------------------------------------
     463             :                         ! Qci_\alpha\alpha'
     464             :                         r_coef_h => jrho1_atom_set(katom)%cjc_ii_h(ispin)%r_coef
     465             :                         r_coef_s => jrho1_atom_set(katom)%cjc_ii_s(ispin)%r_coef
     466             :                         CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
     467             :                                       C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
     468             :                                       c_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
     469             :                                       len_PC1, len_CPC, 1.0_dp, distab)
     470             :                         !------------------------------------------------------------------
     471             :                         ! Qbi_\alpha\alpha'
     472             :                         r_coef_h => jrho1_atom_set(katom)%cjc_iii_h(ispin)%r_coef
     473             :                         r_coef_s => jrho1_atom_set(katom)%cjc_iii_s(ispin)%r_coef
     474             :                         CALL proj_blk(C_coeff_hh_a, C_coeff_ss_a, n_cont_a, &
     475             :                                       C_coeff_hh_b, C_coeff_ss_b, n_cont_b, &
     476             :                                       d_matrix, max_nsgf, r_coef_h, r_coef_s, nso, &
     477             :                                       len_PC1, len_CPC, 1.0_dp, distab)
     478             :                         !------------------------------------------------------------------
     479             : !$                      CALL omp_unset_lock(proj_blk_lock((katom - 1)*nspins + ispin))
     480             : 
     481             :                      END DO ! ispin
     482             : 
     483             :                      EXIT !search loop over jatom-katom list
     484             :                   END IF
     485             :                END DO
     486             :             END DO
     487             : 
     488             :          END DO ! kkind
     489             :          DO ispin = 1, nspins
     490             :             DEALLOCATE (jp_RARnu(ispin)%r_coef, jp2_RARnu(ispin)%r_coef)
     491             :          END DO
     492             :       END DO
     493             :       ! Wait for all threads to finish the loop before locks can be freed
     494             : !$OMP BARRIER
     495             : 
     496             : !$OMP DO
     497             : !$    DO lock = 1, natom
     498             : !$       call omp_destroy_lock(alloc_lock(lock))
     499             : !$    END DO
     500             : !$OMP END DO
     501             : 
     502             : !$OMP DO
     503             : !$    DO lock = 1, nspins*natom
     504             : !$       call omp_destroy_lock(proj_blk_lock(lock))
     505             : !$    END DO
     506             : !$OMP END DO
     507             : 
     508             : !$OMP SINGLE
     509             : !$    DEALLOCATE (alloc_lock)
     510             : !$    DEALLOCATE (proj_blk_lock)
     511             : !$OMP END SINGLE NOWAIT
     512             : 
     513             :       DEALLOCATE (a_matrix, b_matrix, c_matrix, d_matrix, &
     514             :                   a_block, b_block, c_block, d_block, &
     515             :                   jp_RARnu, jp2_RARnu &
     516             :                   )
     517             : 
     518             : !$OMP END PARALLEL
     519             : 
     520         864 :       CALL neighbor_list_iterator_release(nl_iterator)
     521         864 :       DEALLOCATE (basis_set_list)
     522             : 
     523             :       ! parallel sum up
     524         864 :       nbr_dbl = 0.0_dp
     525        2466 :       DO ikind = 1, nkind
     526             :          CALL get_atomic_kind(atomic_kind_set(ikind), &
     527             :                               atom_list=atom_list, &
     528        1602 :                               natom=nat)
     529             :          CALL get_qs_kind(qs_kind_set(ikind), &
     530             :                           basis_set=basis_1c_set, basis_type="GAPW_1C", &
     531        1602 :                           paw_atom=paw_atom)
     532             : 
     533        1602 :          IF (.NOT. paw_atom) CYCLE
     534             : 
     535        1242 :          CALL get_paw_basis_info(basis_1c_set, nsatbas=nsatbas)
     536        1242 :          nsoctot = nsatbas
     537             : 
     538        1242 :          num_pe = para_env%num_pe
     539        1242 :          mepos = para_env%mepos
     540        1242 :          bo = get_limit(nat, num_pe, mepos)
     541             : 
     542        4968 :          ALLOCATE (zero_coeff(nsoctot, nsoctot))
     543        2934 :          DO iat = 1, nat
     544        1692 :             iatom = atom_list(iat)
     545        1692 :             is_not_associated = .NOT. ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(1)%r_coef)
     546             : 
     547        1692 :             IF (iat .GE. bo(1) .AND. iat .LE. bo(2) .AND. is_not_associated) THEN
     548           0 :                CALL allocate_jrho_coeff(jrho1_atom_set, iatom, nsoctot)
     549             :             END IF
     550             : 
     551        5580 :             DO ispin = 1, nspins
     552             : 
     553        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
     554        2646 :                IF (is_not_associated) THEN
     555        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     556             :                END IF
     557     1634454 :                CALL para_env%sum(tmp_coeff)
     558             : 
     559        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
     560        2646 :                IF (is_not_associated) THEN
     561        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     562             :                END IF
     563     1634454 :                CALL para_env%sum(tmp_coeff)
     564             : 
     565        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
     566        2646 :                IF (is_not_associated) THEN
     567        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     568             :                END IF
     569             : 
     570     1634454 :                CALL para_env%sum(tmp_coeff)
     571        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
     572        2646 :                IF (is_not_associated) THEN
     573        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     574             :                END IF
     575     1634454 :                CALL para_env%sum(tmp_coeff)
     576             : 
     577        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
     578        2646 :                IF (is_not_associated) THEN
     579        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     580             :                END IF
     581     1634454 :                CALL para_env%sum(tmp_coeff)
     582             : 
     583        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
     584        2646 :                IF (is_not_associated) THEN
     585        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     586             :                END IF
     587     1634454 :                CALL para_env%sum(tmp_coeff)
     588             : 
     589        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
     590        2646 :                IF (is_not_associated) THEN
     591        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     592             :                END IF
     593     1634454 :                CALL para_env%sum(tmp_coeff)
     594             : 
     595        2646 :                tmp_coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
     596        2646 :                IF (is_not_associated) THEN
     597        1026 :                   zero_coeff = 0.0_dp; tmp_coeff => zero_coeff
     598             :                END IF
     599     1634454 :                CALL para_env%sum(tmp_coeff)
     600             : 
     601        2646 :                IF (ASSOCIATED(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef)) &
     602        4320 :                   nbr_dbl = nbr_dbl + 8.0_dp*REAL(SIZE(jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef), dp)
     603             :             END DO ! ispin
     604             :          END DO ! iat
     605             : 
     606        5310 :          DEALLOCATE (zero_coeff)
     607             : 
     608             :       END DO ! ikind
     609             : 
     610         864 :       output_unit = cp_logger_get_default_io_unit()
     611         864 :       IF (output_unit > 0) THEN
     612         432 :          WRITE (output_unit, '(A,E8.2)') 'calculate_jrho_atom_coeff: nbr_dbl=', nbr_dbl
     613             :       END IF
     614             : 
     615         864 :       CALL timestop(handle)
     616             : 
     617        3456 :    END SUBROUTINE calculate_jrho_atom_coeff
     618             : 
     619             : ! **************************************************************************************************
     620             : !> \brief ...
     621             : !> \param qs_env ...
     622             : !> \param current_env ...
     623             : !> \param idir ...
     624             : ! **************************************************************************************************
     625         864 :    SUBROUTINE calculate_jrho_atom_rad(qs_env, current_env, idir)
     626             :       !
     627             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     628             :       TYPE(current_env_type)                             :: current_env
     629             :       INTEGER, INTENT(IN)                                :: idir
     630             : 
     631             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_jrho_atom_rad'
     632             : 
     633             :       INTEGER :: damax_iso_not0, damax_iso_not0_local, handle, i1, i2, iat, iatom, icg, ikind, &
     634             :          ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso1_first, iso1_last, iso2, iso2_first, &
     635             :          iso2_last, ispin, l, l_iso, llmax, lmax12, lmax_expansion, lmin12, m1s, m2s, m_iso, &
     636             :          max_iso_not0, max_iso_not0_local, max_max_iso_not0, max_nso, max_s_harm, maxl, maxlgto, &
     637             :          maxso, mepos, n1s, n2s, na, natom, natom_tot, nkind, nr, nset, nspins, num_pe, size1, &
     638             :          size2
     639         864 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: cg_n_list, dacg_n_list
     640         864 :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: cg_list, dacg_list
     641             :       INTEGER, DIMENSION(2)                              :: bo
     642         864 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list, lmax, lmin, npgf, o2nindex
     643             :       LOGICAL                                            :: paw_atom
     644         864 :       LOGICAL, ALLOCATABLE, DIMENSION(:, :)              :: is_set_to_0
     645             :       REAL(dp)                                           :: hard_radius
     646         864 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: g1, g2, gauge_h, gauge_s
     647         864 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: cjc0_h_block, cjc0_s_block, cjc_h_block, &
     648         864 :          cjc_ii_h_block, cjc_ii_s_block, cjc_iii_h_block, cjc_iii_s_block, cjc_s_block, dgg_1, gg, &
     649         864 :          gg_lm1
     650         864 :       REAL(dp), DIMENSION(:, :), POINTER :: coeff, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, Fr_a_s, &
     651         864 :          Fr_a_s_ii, Fr_a_s_iii, Fr_b_h, Fr_b_h_ii, Fr_b_h_iii, Fr_b_s, Fr_b_s_ii, Fr_b_s_iii, &
     652         864 :          Fr_h, Fr_s, zet
     653         864 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: my_CG
     654         864 :       REAL(dp), DIMENSION(:, :, :, :), POINTER           :: my_CG_dxyz_asym
     655         864 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     656             :       TYPE(dft_control_type), POINTER                    :: dft_control
     657             :       TYPE(grid_atom_type), POINTER                      :: grid_atom
     658             :       TYPE(gto_basis_set_type), POINTER                  :: basis_1c_set
     659             :       TYPE(harmonics_atom_type), POINTER                 :: harmonics
     660         864 :       TYPE(jrho_atom_type), DIMENSION(:), POINTER        :: jrho1_atom_set
     661             :       TYPE(jrho_atom_type), POINTER                      :: jrho1_atom
     662             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     663         864 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     664             : 
     665         864 :       CALL timeset(routineN, handle)
     666             :       !
     667         864 :       NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, &
     668         864 :                coeff, Fr_h, Fr_s, Fr_a_h, Fr_a_s, Fr_a_h_ii, Fr_a_s_ii, &
     669         864 :                Fr_a_h_iii, Fr_a_s_iii, Fr_b_h, Fr_b_s, Fr_b_h_ii, &
     670         864 :                Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, jrho1_atom_set, &
     671         864 :                jrho1_atom)
     672             :       !
     673             :       CALL get_qs_env(qs_env=qs_env, &
     674             :                       atomic_kind_set=atomic_kind_set, &
     675             :                       qs_kind_set=qs_kind_set, &
     676             :                       dft_control=dft_control, &
     677         864 :                       para_env=para_env)
     678             : 
     679         864 :       CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=maxlgto)
     680             : 
     681             :       !
     682             :       CALL get_current_env(current_env=current_env, &
     683         864 :                            jrho1_atom_set=jrho1_atom_set)
     684             :       !
     685             : 
     686         864 :       nkind = SIZE(qs_kind_set)
     687         864 :       nspins = dft_control%nspins
     688             :       !
     689         864 :       natom_tot = SIZE(jrho1_atom_set, 1)
     690        3456 :       ALLOCATE (is_set_to_0(natom_tot, nspins))
     691        5994 :       is_set_to_0(:, :) = .FALSE.
     692             : 
     693             :       !
     694        2466 :       DO ikind = 1, nkind
     695        1602 :          NULLIFY (atom_list, grid_atom, harmonics, basis_1c_set, &
     696        1602 :                   lmax, lmin, npgf, zet, grid_atom, harmonics, my_CG, my_CG_dxyz_asym)
     697             :          !
     698             :          CALL get_atomic_kind(atomic_kind_set(ikind), &
     699             :                               atom_list=atom_list, &
     700        1602 :                               natom=natom)
     701             :          CALL get_qs_kind(qs_kind_set(ikind), &
     702             :                           grid_atom=grid_atom, &
     703             :                           paw_atom=paw_atom, &
     704             :                           harmonics=harmonics, &
     705             :                           hard_radius=hard_radius, &
     706             :                           basis_set=basis_1c_set, &
     707        1602 :                           basis_type="GAPW_1C")
     708             :          !
     709             :          ! Quick cycle if needed.
     710        1602 :          IF (.NOT. paw_atom) CYCLE
     711             :          !
     712             :          CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
     713             :                                 lmax=lmax, lmin=lmin, &
     714             :                                 maxl=maxl, npgf=npgf, &
     715             :                                 nset=nset, zet=zet, &
     716        1242 :                                 maxso=maxso)
     717        1242 :          CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex)
     718             :          !
     719        1242 :          nr = grid_atom%nr
     720        1242 :          na = grid_atom%ng_sphere
     721        1242 :          max_iso_not0 = harmonics%max_iso_not0
     722        1242 :          damax_iso_not0 = harmonics%damax_iso_not0
     723        1242 :          max_max_iso_not0 = MAX(max_iso_not0, damax_iso_not0)
     724        1242 :          lmax_expansion = indso(1, max_max_iso_not0)
     725        1242 :          max_s_harm = harmonics%max_s_harm
     726        1242 :          llmax = harmonics%llmax
     727             :          !
     728             :          ! Distribute the atoms of this kind
     729        1242 :          num_pe = para_env%num_pe
     730        1242 :          mepos = para_env%mepos
     731        1242 :          bo = get_limit(natom, num_pe, mepos)
     732             :          !
     733        1242 :          my_CG => harmonics%my_CG
     734        1242 :          my_CG_dxyz_asym => harmonics%my_CG_dxyz_asym
     735             :          !
     736             :          ! Allocate some arrays.
     737        1242 :          max_nso = nsoset(maxl)
     738           0 :          ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), gg_lm1(nr, 0:2*maxl), dgg_1(nr, 0:2*maxl), &
     739           0 :                    cjc0_h_block(max_nso, max_nso), cjc0_s_block(max_nso, max_nso), &
     740           0 :                    cjc_h_block(max_nso, max_nso), cjc_s_block(max_nso, max_nso), &
     741           0 :                    cjc_ii_h_block(max_nso, max_nso), cjc_ii_s_block(max_nso, max_nso), &
     742           0 :                    cjc_iii_h_block(max_nso, max_nso), cjc_iii_s_block(max_nso, max_nso), &
     743           0 :                    cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
     744           0 :                    dacg_list(2, nsoset(maxl)**2, max_s_harm), dacg_n_list(max_s_harm), &
     745       64584 :                    gauge_h(nr), gauge_s(nr))
     746             :          !
     747             :          ! Compute the gauge
     748        1296 :          SELECT CASE (current_env%gauge)
     749             :          CASE (current_gauge_r)
     750             :             ! d(r)=r
     751        2754 :             gauge_h(1:nr) = grid_atom%rad(1:nr)
     752        2754 :             gauge_s(1:nr) = grid_atom%rad(1:nr)
     753             :          CASE (current_gauge_r_and_step_func)
     754             :             ! step function
     755       43740 :             gauge_h(1:nr) = 0e0_dp
     756       43740 :             DO ir = 1, nr
     757       43740 :                IF (grid_atom%rad(ir) .LE. hard_radius) THEN
     758       24138 :                   gauge_s(ir) = grid_atom%rad(ir)
     759             :                ELSE
     760       18702 :                   gauge_s(ir) = gauge_h(ir)
     761             :                END IF
     762             :             END DO
     763             :          CASE (current_gauge_atom)
     764             :             ! d(r)=A
     765       15588 :             gauge_h(1:nr) = HUGE(0e0_dp) !0e0_dp
     766       15588 :             gauge_s(1:nr) = HUGE(0e0_dp) !0e0_dp
     767             :          CASE DEFAULT
     768        1242 :             CPABORT("Unknown gauge, try again...")
     769             :          END SELECT
     770             :          !
     771             :          !
     772        1242 :          m1s = 0
     773        4968 :          DO iset1 = 1, nset
     774             :             m2s = 0
     775       15300 :             DO iset2 = 1, nset
     776             :                CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
     777       11574 :                                       max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
     778       11574 :                CPASSERT(max_iso_not0_local .LE. max_iso_not0)
     779             :                CALL get_none0_cg_list(my_CG_dxyz_asym, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
     780       11574 :                                       max_s_harm, lmax_expansion, dacg_list, dacg_n_list, damax_iso_not0_local)
     781       11574 :                CPASSERT(damax_iso_not0_local .LE. damax_iso_not0)
     782             : 
     783       11574 :                n1s = nsoset(lmax(iset1))
     784       36504 :                DO ipgf1 = 1, npgf(iset1)
     785       24930 :                   iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
     786       24930 :                   iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
     787       24930 :                   size1 = iso1_last - iso1_first + 1
     788       24930 :                   iso1_first = o2nindex(iso1_first)
     789       24930 :                   iso1_last = o2nindex(iso1_last)
     790       24930 :                   i1 = iso1_last - iso1_first + 1
     791       24930 :                   CPASSERT(size1 == i1)
     792       24930 :                   i1 = nsoset(lmin(iset1) - 1) + 1
     793             :                   !
     794     1244430 :                   g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
     795             :                   !
     796       24930 :                   n2s = nsoset(lmax(iset2))
     797       91674 :                   DO ipgf2 = 1, npgf(iset2)
     798       55170 :                      iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
     799       55170 :                      iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
     800       55170 :                      size2 = iso2_last - iso2_first + 1
     801       55170 :                      iso2_first = o2nindex(iso2_first)
     802       55170 :                      iso2_last = o2nindex(iso2_last)
     803       55170 :                      i2 = iso2_last - iso2_first + 1
     804       55170 :                      CPASSERT(size2 == i2)
     805       55170 :                      i2 = nsoset(lmin(iset2) - 1) + 1
     806             :                      !
     807     2730510 :                      g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
     808             :                      !
     809       55170 :                      lmin12 = lmin(iset1) + lmin(iset2)
     810       55170 :                      lmax12 = lmax(iset1) + lmax(iset2)
     811             :                      !
     812    10133028 :                      gg = 0.0_dp
     813    10133028 :                      gg_lm1 = 0.0_dp
     814    10133028 :                      dgg_1 = 0.0_dp
     815             :                      !
     816             :                      ! Take only the terms of expansion for L < lmax_expansion
     817       55170 :                      IF (lmin12 .LE. lmax_expansion) THEN
     818             :                         !
     819             :                         IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
     820             :                         !
     821       55170 :                         IF (lmin12 == 0) THEN
     822     2510514 :                            gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
     823     2510514 :                            gg_lm1(1:nr, lmin12) = 0.0_dp
     824             :                         ELSE
     825      219996 :                            gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
     826      219996 :                            gg_lm1(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12 - 1)*g1(1:nr)*g2(1:nr)
     827             :                         END IF
     828             :                         !
     829      103122 :                         DO l = lmin12 + 1, lmax12
     830     2363472 :                            gg(1:nr, l) = grid_atom%rad(1:nr)*gg(1:nr, l - 1)
     831     2418642 :                            gg_lm1(1:nr, l) = gg(1:nr, l - 1)
     832             :                         END DO
     833             :                         !
     834      158292 :                         DO l = lmin12, lmax12
     835             :                            dgg_1(1:nr, l) = 2.0_dp*(zet(ipgf1, iset1) - zet(ipgf2, iset2))&
     836     5149152 :                                            &              *gg(1:nr, l)*grid_atom%rad(1:nr)
     837             :                         END DO
     838             :                      ELSE
     839             :                         CYCLE
     840             :                      END IF ! lmin12
     841             :                      !
     842      118251 :                      DO iat = bo(1), bo(2)
     843       38151 :                         iatom = atom_list(iat)
     844             :                         !
     845      155052 :                         DO ispin = 1, nspins
     846             :                            !------------------------------------------------------------------
     847             :                            ! P_\alpha\alpha'
     848     3125673 :                            cjc0_h_block = HUGE(1.0_dp)
     849     3125673 :                            cjc0_s_block = HUGE(1.0_dp)
     850             :                            !
     851             :                            ! Hard term
     852       61731 :                            coeff => jrho1_atom_set(iatom)%cjc0_h(ispin)%r_coef
     853             :                            cjc0_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     854      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     855             :                            !
     856             :                            ! Soft term
     857       61731 :                            coeff => jrho1_atom_set(iatom)%cjc0_s(ispin)%r_coef
     858             :                            cjc0_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     859      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     860             :                            !------------------------------------------------------------------
     861             :                            ! mQai_\alpha\alpha'
     862     3125673 :                            cjc_h_block = HUGE(1.0_dp)
     863     3125673 :                            cjc_s_block = HUGE(1.0_dp)
     864             :                            !
     865             :                            ! Hard term
     866       61731 :                            coeff => jrho1_atom_set(iatom)%cjc_h(ispin)%r_coef
     867             :                            cjc_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     868      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     869             :                            !
     870             :                            ! Soft term
     871       61731 :                            coeff => jrho1_atom_set(iatom)%cjc_s(ispin)%r_coef
     872             :                            cjc_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     873      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     874             :                            !------------------------------------------------------------------
     875             :                            ! Qci_\alpha\alpha'
     876     3125673 :                            cjc_ii_h_block = HUGE(1.0_dp)
     877     3125673 :                            cjc_ii_s_block = HUGE(1.0_dp)
     878             :                            !
     879             :                            ! Hard term
     880       61731 :                            coeff => jrho1_atom_set(iatom)%cjc_ii_h(ispin)%r_coef
     881             :                            cjc_ii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     882      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     883             :                            !
     884             :                            ! Soft term
     885       61731 :                            coeff => jrho1_atom_set(iatom)%cjc_ii_s(ispin)%r_coef
     886             :                            cjc_ii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     887      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     888             :                            !------------------------------------------------------------------
     889             :                            ! Qbi_\alpha\alpha'
     890     3125673 :                            cjc_iii_h_block = HUGE(1.0_dp)
     891     3125673 :                            cjc_iii_s_block = HUGE(1.0_dp)
     892             :                            !
     893             :                            !
     894             :                            ! Hard term
     895       61731 :                            coeff => jrho1_atom_set(iatom)%cjc_iii_h(ispin)%r_coef
     896             :                            cjc_iii_h_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     897      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     898             :                            !
     899             :                            ! Soft term
     900       61731 :                            coeff => jrho1_atom_set(iatom)%cjc_iii_s(ispin)%r_coef
     901             :                            cjc_iii_s_block(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
     902      601002 :                               coeff(iso1_first:iso1_last, iso2_first:iso2_last)
     903             :                            !------------------------------------------------------------------
     904             :                            !
     905             :                            ! Allocation radial functions
     906       61731 :                            jrho1_atom => jrho1_atom_set(iatom)
     907       61731 :                            IF (.NOT. ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef)) THEN
     908             :                               CALL allocate_jrho_atom_rad(jrho1_atom, ispin, nr, na, &
     909         147 :                                                           max_max_iso_not0)
     910         147 :                               is_set_to_0(iatom, ispin) = .TRUE.
     911             :                            ELSE
     912       61584 :                               IF (.NOT. is_set_to_0(iatom, ispin)) THEN
     913        1176 :                                  CALL set2zero_jrho_atom_rad(jrho1_atom, ispin)
     914        1176 :                                  is_set_to_0(iatom, ispin) = .TRUE.
     915             :                               END IF
     916             :                            END IF
     917             :                            !------------------------------------------------------------------
     918             :                            !
     919       61731 :                            Fr_h => jrho1_atom%jrho_h(ispin)%r_coef
     920       61731 :                            Fr_s => jrho1_atom%jrho_s(ispin)%r_coef
     921             :                            !------------------------------------------------------------------
     922             :                            !
     923       61731 :                            Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef
     924       61731 :                            Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
     925       61731 :                            Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef
     926       61731 :                            Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
     927             :                            !------------------------------------------------------------------
     928             :                            !
     929       61731 :                            Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef
     930       61731 :                            Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
     931       61731 :                            Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef
     932       61731 :                            Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
     933             :                            !------------------------------------------------------------------
     934             :                            !
     935       61731 :                            Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef
     936       61731 :                            Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
     937       61731 :                            Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef
     938       61731 :                            Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
     939             :                            !------------------------------------------------------------------
     940             :                            !
     941      362160 :                            DO iso = 1, max_iso_not0_local
     942      300429 :                               l_iso = indso(1, iso) ! not needed
     943      300429 :                               m_iso = indso(2, iso) ! not needed
     944             :                               !
     945      858771 :                               DO icg = 1, cg_n_list(iso)
     946             :                                  !
     947      496611 :                                  iso1 = cg_list(1, icg, iso)
     948      496611 :                                  iso2 = cg_list(2, icg, iso)
     949             :                                  !
     950      496611 :                                  IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
     951           0 :                                     WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
     952           0 :                                     WRITE (*, *) '.... will stop!'
     953             :                                  END IF
     954      496611 :                                  CPASSERT(iso2 > 0 .AND. iso1 > 0)
     955             :                                  !
     956      496611 :                                  l = indso(1, iso1) + indso(1, iso2)
     957      496611 :                                  IF (l .GT. lmax_expansion .OR. l .LT. .0) THEN
     958           0 :                                     WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
     959           0 :                                     WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
     960           0 :                                     WRITE (*, *) '.... will stop!'
     961             :                                  END IF
     962      496611 :                                  CPASSERT(l <= lmax_expansion)
     963             :                                  !------------------------------------------------------------------
     964             :                                  ! P0
     965             :                                  !
     966      496611 :                                  IF (current_env%gauge .EQ. current_gauge_atom) THEN
     967             :                                     ! Hard term
     968             :                                     Fr_h(1:nr, iso) = Fr_h(1:nr, iso) + &
     969             :                                                       gg(1:nr, l)*cjc0_h_block(iso1, iso2)* &
     970     4477842 :                                                       my_CG(iso1, iso2, iso)
     971             :                                     ! Soft term
     972             :                                     Fr_s(1:nr, iso) = Fr_s(1:nr, iso) + &
     973             :                                                       gg(1:nr, l)*cjc0_s_block(iso1, iso2)* &
     974     4477842 :                                                       my_CG(iso1, iso2, iso)
     975             :                                  ELSE
     976             :                                     ! Hard term
     977             :                                     Fr_h(1:nr, iso) = Fr_h(1:nr, iso) + &
     978             :                                                       gg(1:nr, l)*cjc0_h_block(iso1, iso2)* &
     979    39661029 :                                                       my_CG(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_h(1:nr))
     980             :                                     ! Soft term
     981             :                                     Fr_s(1:nr, iso) = Fr_s(1:nr, iso) + &
     982             :                                                       gg(1:nr, l)*cjc0_s_block(iso1, iso2)* &
     983    39661029 :                                                       my_CG(iso1, iso2, iso)*(grid_atom%rad(1:nr) - gauge_s(1:nr))
     984             :                                  END IF
     985             :                                  !------------------------------------------------------------------
     986             :                                  ! Rai
     987             :                                  !
     988             :                                  ! Hard term
     989             :                                  Fr_a_h(1:nr, iso) = Fr_a_h(1:nr, iso) + &
     990             :                                                      dgg_1(1:nr, l)*cjc_h_block(iso1, iso2)* &
     991    24512841 :                                                      my_CG(iso1, iso2, iso)
     992             :                                  !
     993             :                                  ! Soft term
     994             :                                  Fr_a_s(1:nr, iso) = Fr_a_s(1:nr, iso) + &
     995             :                                                      dgg_1(1:nr, l)*cjc_s_block(iso1, iso2)* &
     996    24512841 :                                                      my_CG(iso1, iso2, iso)
     997             :                                  !------------------------------------------------------------------
     998             :                                  ! Rci
     999             :                                  !
    1000      496611 :                                  IF (current_env%gauge .EQ. current_gauge_atom) THEN
    1001             :                                     ! Hard term
    1002             :                                     Fr_a_h_ii(1:nr, iso) = Fr_a_h_ii(1:nr, iso) + &
    1003             :                                                            dgg_1(1:nr, l)* &
    1004             :                                                            cjc_ii_h_block(iso1, iso2)* &
    1005     4477842 :                                                            my_CG(iso1, iso2, iso)
    1006             :                                     ! Soft term
    1007             :                                     Fr_a_s_ii(1:nr, iso) = Fr_a_s_ii(1:nr, iso) + &
    1008             :                                                            dgg_1(1:nr, l)* &
    1009             :                                                            cjc_ii_s_block(iso1, iso2)* &
    1010     4477842 :                                                            my_CG(iso1, iso2, iso)
    1011             :                                  ELSE
    1012             :                                     ! Hard term
    1013             :                                     Fr_a_h_ii(1:nr, iso) = Fr_a_h_ii(1:nr, iso) + &
    1014             :                                                            dgg_1(1:nr, l)*gauge_h(1:nr)* &
    1015             :                                                            cjc_ii_h_block(iso1, iso2)* &
    1016    20034999 :                                                            my_CG(iso1, iso2, iso)
    1017             :                                     ! Soft term
    1018             :                                     Fr_a_s_ii(1:nr, iso) = Fr_a_s_ii(1:nr, iso) + &
    1019             :                                                            dgg_1(1:nr, l)*gauge_s(1:nr)* &
    1020             :                                                            cjc_ii_s_block(iso1, iso2)* &
    1021    20034999 :                                                            my_CG(iso1, iso2, iso)
    1022             :                                  END IF
    1023             :                                  !------------------------------------------------------------------
    1024             :                                  ! Rbi
    1025             :                                  !
    1026      797040 :                                  IF (current_env%gauge .EQ. current_gauge_atom) THEN
    1027             :                                     ! Hard term
    1028             :                                     Fr_a_h_iii(1:nr, iso) = Fr_a_h_iii(1:nr, iso) + &
    1029             :                                                             dgg_1(1:nr, l)* &
    1030             :                                                             cjc_iii_h_block(iso1, iso2)* &
    1031     4477842 :                                                             my_CG(iso1, iso2, iso)
    1032             :                                     ! Soft term
    1033             :                                     Fr_a_s_iii(1:nr, iso) = Fr_a_s_iii(1:nr, iso) + &
    1034             :                                                             dgg_1(1:nr, l)* &
    1035             :                                                             cjc_iii_s_block(iso1, iso2)* &
    1036     4477842 :                                                             my_CG(iso1, iso2, iso)
    1037             :                                  ELSE
    1038             :                                     ! Hard term
    1039             :                                     Fr_a_h_iii(1:nr, iso) = Fr_a_h_iii(1:nr, iso) + &
    1040             :                                                             dgg_1(1:nr, l)*gauge_h(1:nr)* &
    1041             :                                                             cjc_iii_h_block(iso1, iso2)* &
    1042    20034999 :                                                             my_CG(iso1, iso2, iso)
    1043             :                                     ! Soft term
    1044             :                                     Fr_a_s_iii(1:nr, iso) = Fr_a_s_iii(1:nr, iso) + &
    1045             :                                                             dgg_1(1:nr, l)*gauge_s(1:nr)* &
    1046             :                                                             cjc_iii_s_block(iso1, iso2)* &
    1047    20034999 :                                                             my_CG(iso1, iso2, iso)
    1048             :                                  END IF
    1049             :                                  !------------------------------------------------------------------
    1050             :                               END DO !icg
    1051             :                               !
    1052             :                            END DO ! iso
    1053             :                            !
    1054             :                            !
    1055      212436 :                            DO iso = 1, damax_iso_not0_local
    1056      112554 :                               l_iso = indso(1, iso) ! not needed
    1057      112554 :                               m_iso = indso(2, iso) ! not needed
    1058             :                               !
    1059      669321 :                               DO icg = 1, dacg_n_list(iso)
    1060             :                                  !
    1061      495036 :                                  iso1 = dacg_list(1, icg, iso)
    1062      495036 :                                  iso2 = dacg_list(2, icg, iso)
    1063             :                                  !
    1064      495036 :                                  IF (.NOT. (iso2 > 0 .AND. iso1 > 0)) THEN
    1065           0 :                                     WRITE (*, *) 'iso1=', iso1, ' iso2=', iso2, ' iso=', iso, ' icg=', icg
    1066           0 :                                     WRITE (*, *) '.... will stop!'
    1067             :                                  END IF
    1068      495036 :                                  CPASSERT(iso2 > 0 .AND. iso1 > 0)
    1069             :                                  !
    1070      495036 :                                  l = indso(1, iso1) + indso(1, iso2)
    1071      495036 :                                  IF (l .GT. lmax_expansion) THEN
    1072           0 :                                     WRITE (*, *) 'calculate_jrho_atom_rad: 1 l', l
    1073           0 :                                     WRITE (*, *) 'calculate_jrho_atom_rad: 1 lmax_expansion', lmax_expansion
    1074           0 :                                     WRITE (*, *) '.... will stop!'
    1075             :                                  END IF
    1076      495036 :                                  CPASSERT(l <= lmax_expansion)
    1077             :                                  !------------------------------------------------------------------
    1078             :                                  ! Daij
    1079             :                                  !
    1080             :                                  ! Hard term
    1081             :                                  Fr_b_h(1:nr, iso) = Fr_b_h(1:nr, iso) + &
    1082             :                                                      gg_lm1(1:nr, l)*cjc_h_block(iso1, iso2)* &
    1083    24139836 :                                                      my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1084             :                                  !
    1085             :                                  ! Soft term
    1086             :                                  Fr_b_s(1:nr, iso) = Fr_b_s(1:nr, iso) + &
    1087             :                                                      gg_lm1(1:nr, l)*cjc_s_block(iso1, iso2)* &
    1088    24139836 :                                                      my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1089             :                                  !
    1090             :                                  !------------------------------------------------------------------
    1091             :                                  ! Dcij
    1092             :                                  !
    1093      495036 :                                  IF (current_env%gauge .EQ. current_gauge_atom) THEN
    1094             :                                     ! Hard term
    1095             :                                     Fr_b_h_ii(1:nr, iso) = Fr_b_h_ii(1:nr, iso) + &
    1096             :                                                            gg_lm1(1:nr, l)* &
    1097             :                                                            cjc_ii_h_block(iso1, iso2)* &
    1098     3569184 :                                                            my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1099             :                                     ! Soft term
    1100             :                                     Fr_b_s_ii(1:nr, iso) = Fr_b_s_ii(1:nr, iso) + &
    1101             :                                                            gg_lm1(1:nr, l)* &
    1102             :                                                            cjc_ii_s_block(iso1, iso2)* &
    1103     3569184 :                                                            my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1104             :                                  ELSE
    1105             :                                     ! Hard term
    1106             :                                     Fr_b_h_ii(1:nr, iso) = Fr_b_h_ii(1:nr, iso) + &
    1107             :                                                            gg_lm1(1:nr, l)*gauge_h(1:nr)* &
    1108             :                                                            cjc_ii_h_block(iso1, iso2)* &
    1109    20570652 :                                                            my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1110             :                                     ! Soft term
    1111             :                                     Fr_b_s_ii(1:nr, iso) = Fr_b_s_ii(1:nr, iso) + &
    1112             :                                                            gg_lm1(1:nr, l)*gauge_s(1:nr)* &
    1113             :                                                            cjc_ii_s_block(iso1, iso2)* &
    1114    20570652 :                                                            my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1115             :                                  END IF
    1116             :                                  !------------------------------------------------------------------
    1117             :                                  ! Dbij
    1118             :                                  !
    1119      607590 :                                  IF (current_env%gauge .EQ. current_gauge_atom) THEN
    1120             :                                     ! Hard term
    1121             :                                     Fr_b_h_iii(1:nr, iso) = Fr_b_h_iii(1:nr, iso) + &
    1122             :                                                             gg_lm1(1:nr, l)* &
    1123             :                                                             cjc_iii_h_block(iso1, iso2)* &
    1124     3569184 :                                                             my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1125             :                                     ! Soft term
    1126             :                                     Fr_b_s_iii(1:nr, iso) = Fr_b_s_iii(1:nr, iso) + &
    1127             :                                                             gg_lm1(1:nr, l)* &
    1128             :                                                             cjc_iii_s_block(iso1, iso2)* &
    1129     3569184 :                                                             my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1130             :                                  ELSE
    1131             :                                     ! Hard term
    1132             :                                     Fr_b_h_iii(1:nr, iso) = Fr_b_h_iii(1:nr, iso) + &
    1133             :                                                             gg_lm1(1:nr, l)*gauge_h(1:nr)* &
    1134             :                                                             cjc_iii_h_block(iso1, iso2)* &
    1135    20570652 :                                                             my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1136             :                                     ! Soft term
    1137             :                                     Fr_b_s_iii(1:nr, iso) = Fr_b_s_iii(1:nr, iso) + &
    1138             :                                                             gg_lm1(1:nr, l)*gauge_s(1:nr)* &
    1139             :                                                             cjc_iii_s_block(iso1, iso2)* &
    1140    20570652 :                                                             my_CG_dxyz_asym(idir, iso1, iso2, iso)
    1141             :                                  END IF
    1142             :                                  !------------------------------------------------------------------
    1143             :                               END DO ! icg
    1144             :                            END DO ! iso
    1145             :                            !
    1146             :                         END DO ! ispin
    1147             :                      END DO ! iat
    1148             :                      !
    1149             :                      !------------------------------------------------------------------
    1150             :                      !
    1151             :                   END DO !ipgf2
    1152             :                END DO ! ipgf1
    1153       38448 :                m2s = m2s + maxso
    1154             :             END DO ! iset2
    1155        4968 :             m1s = m1s + maxso
    1156             :          END DO ! iset1
    1157             :          !
    1158           0 :          DEALLOCATE (cjc0_h_block, cjc0_s_block, cjc_h_block, cjc_s_block, cjc_ii_h_block, cjc_ii_s_block, &
    1159           0 :                      cjc_iii_h_block, cjc_iii_s_block, g1, g2, gg, gg_lm1, dgg_1, gauge_h, gauge_s, &
    1160        1242 :                      cg_list, cg_n_list, dacg_list, dacg_n_list)
    1161        5310 :          DEALLOCATE (o2nindex)
    1162             :       END DO ! ikind
    1163             :       !
    1164             :       !
    1165         864 :       DEALLOCATE (is_set_to_0)
    1166             :       !
    1167         864 :       CALL timestop(handle)
    1168             :       !
    1169        2592 :    END SUBROUTINE calculate_jrho_atom_rad
    1170             : 
    1171             : ! **************************************************************************************************
    1172             : !> \brief ...
    1173             : !> \param jrho1_atom ...
    1174             : !> \param jrho_h ...
    1175             : !> \param jrho_s ...
    1176             : !> \param grid_atom ...
    1177             : !> \param harmonics ...
    1178             : !> \param do_igaim ...
    1179             : !> \param ratom ...
    1180             : !> \param natm_gauge ...
    1181             : !> \param iB ...
    1182             : !> \param idir ...
    1183             : !> \param ispin ...
    1184             : ! **************************************************************************************************
    1185        1323 :    SUBROUTINE calculate_jrho_atom_ang(jrho1_atom, jrho_h, jrho_s, grid_atom, &
    1186        1323 :                                       harmonics, do_igaim, ratom, natm_gauge, &
    1187             :                                       iB, idir, ispin)
    1188             :       !
    1189             :       TYPE(jrho_atom_type), POINTER            :: jrho1_atom
    1190             :       REAL(dp), DIMENSION(:, :), POINTER       :: jrho_h, jrho_s
    1191             :       TYPE(grid_atom_type), POINTER            :: grid_atom
    1192             :       TYPE(harmonics_atom_type), POINTER       :: harmonics
    1193             :       LOGICAL, INTENT(IN)                      :: do_igaim
    1194             :       INTEGER, INTENT(IN)                      :: iB, idir, ispin, natm_gauge
    1195             :       REAL(dp), INTENT(IN) :: ratom(:, :)
    1196             : 
    1197             :       INTEGER                                  :: ia, idir2, iiB, iiiB, ir, &
    1198             :                                                   iso, max_max_iso_not0, na, nr
    1199             :       REAL(dp)                                 :: rad_part, scale
    1200        1323 :       REAL(dp), DIMENSION(:, :), POINTER :: a, Fr_a_h, Fr_a_h_ii, Fr_a_h_iii, &
    1201        1323 :                                             Fr_a_s, Fr_a_s_ii, Fr_a_s_iii, Fr_b_h, Fr_b_h_ii, Fr_b_h_iii, Fr_b_s, &
    1202        1323 :                                             Fr_b_s_ii, Fr_b_s_iii, Fr_h, Fr_s, slm
    1203        1323 :       REAL(dp), DIMENSION(:), POINTER :: r
    1204        1323 :       REAL(dp), DIMENSION(:, :, :), ALLOCATABLE :: g
    1205             : !
    1206             : !
    1207        1323 :       NULLIFY (Fr_h, Fr_s, Fr_a_h, Fr_a_s, Fr_a_h_ii, Fr_a_s_ii, Fr_a_h_iii, Fr_a_s_iii, &
    1208        1323 :                Fr_b_h, Fr_b_s, Fr_b_h_ii, Fr_b_s_ii, Fr_b_h_iii, Fr_b_s_iii, &
    1209        1323 :                a, slm)
    1210             :       !
    1211           0 :       CPASSERT(ASSOCIATED(jrho_h))
    1212        1323 :       CPASSERT(ASSOCIATED(jrho_s))
    1213        1323 :       CPASSERT(ASSOCIATED(jrho1_atom))
    1214             :       ! just to be sure...
    1215        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h))
    1216        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s))
    1217        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h))
    1218        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s))
    1219        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h(ispin)%r_coef))
    1220        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s(ispin)%r_coef))
    1221        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h(ispin)%r_coef))
    1222        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s(ispin)%r_coef))
    1223        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii))
    1224        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii))
    1225        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii))
    1226        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii))
    1227        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_ii(ispin)%r_coef))
    1228        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_ii(ispin)%r_coef))
    1229        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_ii(ispin)%r_coef))
    1230        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_ii(ispin)%r_coef))
    1231        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii))
    1232        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii))
    1233        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii))
    1234        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii))
    1235        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_h_iii(ispin)%r_coef))
    1236        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_a_s_iii(ispin)%r_coef))
    1237        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_h_iii(ispin)%r_coef))
    1238        1323 :       CPASSERT(ASSOCIATED(jrho1_atom%jrho_b_s_iii(ispin)%r_coef))
    1239             :       !
    1240             :       !
    1241        1323 :       nr = grid_atom%nr
    1242        1323 :       na = grid_atom%ng_sphere
    1243        1323 :       max_max_iso_not0 = MAX(harmonics%max_iso_not0, harmonics%damax_iso_not0)
    1244        5292 :       ALLOCATE (g(3, nr, na))
    1245             :       !------------------------------------------------------------------
    1246             :       !
    1247        1323 :       Fr_h => jrho1_atom%jrho_h(ispin)%r_coef
    1248        1323 :       Fr_s => jrho1_atom%jrho_s(ispin)%r_coef
    1249             :       !------------------------------------------------------------------
    1250             :       !
    1251        1323 :       Fr_a_h => jrho1_atom%jrho_a_h(ispin)%r_coef !Rai
    1252        1323 :       Fr_a_s => jrho1_atom%jrho_a_s(ispin)%r_coef
    1253        1323 :       Fr_b_h => jrho1_atom%jrho_b_h(ispin)%r_coef !Daij
    1254        1323 :       Fr_b_s => jrho1_atom%jrho_b_s(ispin)%r_coef
    1255             :       !------------------------------------------------------------------
    1256             :       !
    1257        1323 :       Fr_a_h_ii => jrho1_atom%jrho_a_h_ii(ispin)%r_coef !Rci
    1258        1323 :       Fr_a_s_ii => jrho1_atom%jrho_a_s_ii(ispin)%r_coef
    1259        1323 :       Fr_b_h_ii => jrho1_atom%jrho_b_h_ii(ispin)%r_coef !Dcij
    1260        1323 :       Fr_b_s_ii => jrho1_atom%jrho_b_s_ii(ispin)%r_coef
    1261             :       !------------------------------------------------------------------
    1262             :       !
    1263        1323 :       Fr_a_h_iii => jrho1_atom%jrho_a_h_iii(ispin)%r_coef !Rbi
    1264        1323 :       Fr_a_s_iii => jrho1_atom%jrho_a_s_iii(ispin)%r_coef
    1265        1323 :       Fr_b_h_iii => jrho1_atom%jrho_b_h_iii(ispin)%r_coef !Dbij
    1266        1323 :       Fr_b_s_iii => jrho1_atom%jrho_b_s_iii(ispin)%r_coef
    1267             :       !------------------------------------------------------------------
    1268             :       !
    1269        1323 :       a => harmonics%a
    1270        1323 :       slm => harmonics%slm
    1271        1323 :       r => grid_atom%rad
    1272             :       !
    1273        1323 :       CALL set_vecp(iB, iiB, iiiB)
    1274             :       !
    1275        1323 :       scale = 0.0_dp
    1276        1323 :       idir2 = 1
    1277        1323 :       IF (idir .NE. iB) THEN
    1278         882 :          CALL set_vecp_rev(idir, iB, idir2)
    1279         882 :          scale = fac_vecp(idir, iB, idir2)
    1280             :       END IF
    1281             :       !
    1282             :       ! Set the gauge
    1283        1323 :       CALL get_gauge()
    1284             :       !
    1285       65943 :       DO ir = 1, nr
    1286      820323 :          DO iso = 1, max_max_iso_not0
    1287    38013120 :             DO ia = 1, na
    1288    37948500 :                IF (do_igaim) THEN
    1289             :                   !------------------------------------------------------------------
    1290             :                   ! Hard current density response
    1291             :                   ! radial(ia,ir) = (               aj(ia) * Rai(ir,iso) + Daij
    1292             :                   !                  -  aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij )
    1293             :                   !                  + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij )
    1294             :                   !                 ) * Ylm(ia)
    1295             :                   rad_part = a(idir, ia)*Fr_a_h(ir, iso) + Fr_b_h(ir, iso) &
    1296             :                        &   - g(iiB, ir, ia)*(a(idir, ia)*Fr_a_h_iii(ir, iso) + Fr_b_h_iii(ir, iso))&
    1297             :                        &   + g(iiiB, ir, ia)*(a(idir, ia)*Fr_a_h_ii(ir, iso) + Fr_b_h_ii(ir, iso))&
    1298     7380000 :                        &   + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_h(ir, iso)
    1299             :                   !
    1300     7380000 :                   jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso)
    1301             :                   !------------------------------------------------------------------
    1302             :                   ! Soft current density response
    1303             :                   rad_part = a(idir, ia)*Fr_a_s(ir, iso) + Fr_b_s(ir, iso) &
    1304             :                        &   - g(iiB, ir, ia)*(a(idir, ia)*Fr_a_s_iii(ir, iso) + Fr_b_s_iii(ir, iso))&
    1305             :                        &   + g(iiiB, ir, ia)*(a(idir, ia)*Fr_a_s_ii(ir, iso) + Fr_b_s_ii(ir, iso))&
    1306     7380000 :                        &   + scale*(a(idir2, ia)*r(ir) - g(idir2, ir, ia))*Fr_s(ir, iso)
    1307             :                   !
    1308     7380000 :                   jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso)
    1309             :                   !------------------------------------------------------------------
    1310             :                ELSE
    1311             :                   !------------------------------------------------------------------
    1312             :                   ! Hard current density response
    1313             :                   ! radial(ia,ir) = (               aj(ia) * Rai(ir,iso) + Daij
    1314             :                   !                  -  aii(ia) * ( aj(ia) * Rbi(ir,iso) + Dbij )
    1315             :                   !                  + aiii(ia) * ( aj(ia) * Rci(ir,iso) + Dcij )
    1316             :                   !                 ) * Ylm(ia)
    1317             :                   rad_part = a(idir, ia)*Fr_a_h(ir, iso) + Fr_b_h(ir, iso) &
    1318             :                        &   - a(iiB, ia)*(a(idir, ia)*Fr_a_h_iii(ir, iso) + Fr_b_h_iii(ir, iso))&
    1319             :                        &   + a(iiiB, ia)*(a(idir, ia)*Fr_a_h_ii(ir, iso) + Fr_b_h_ii(ir, iso))&
    1320    29814120 :                        &   + scale*a(idir2, ia)*Fr_h(ir, iso)
    1321             :                   !
    1322    29814120 :                   jrho_h(ir, ia) = jrho_h(ir, ia) + rad_part*slm(ia, iso)
    1323             :                   !------------------------------------------------------------------
    1324             :                   ! Soft current density response
    1325             :                   rad_part = a(idir, ia)*Fr_a_s(ir, iso) + Fr_b_s(ir, iso) &
    1326             :                        &   - a(iiB, ia)*(a(idir, ia)*Fr_a_s_iii(ir, iso) + Fr_b_s_iii(ir, iso))&
    1327             :                        &   + a(iiiB, ia)*(a(idir, ia)*Fr_a_s_ii(ir, iso) + Fr_b_s_ii(ir, iso))&
    1328    29814120 :                        &   + scale*a(idir2, ia)*Fr_s(ir, iso)
    1329             :                   !
    1330    29814120 :                   jrho_s(ir, ia) = jrho_s(ir, ia) + rad_part*slm(ia, iso)
    1331             :                   !------------------------------------------------------------------
    1332             :                END IF
    1333             :             END DO ! ia
    1334             :          END DO ! iso
    1335             :       END DO ! ir
    1336             :       !
    1337        3969 :       DEALLOCATE (g)
    1338             :       !
    1339             :    CONTAINS
    1340             :       !
    1341             : ! **************************************************************************************************
    1342             : !> \brief ...
    1343             : ! **************************************************************************************************
    1344        1323 :       SUBROUTINE get_gauge()
    1345             :       INTEGER                                            :: iatom, ixyz, jatom
    1346             :       REAL(dp)                                           :: ab, pa, pb, point(3), pra(3), prb(3), &
    1347             :                                                             res, tmp
    1348        1323 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: buf
    1349             : 
    1350        3969 :          ALLOCATE (buf(natm_gauge))
    1351       65943 :          DO ir = 1, nr
    1352     3238623 :          DO ia = 1, na
    1353    12690720 :             DO ixyz = 1, 3
    1354    12690720 :                g(ixyz, ir, ia) = 0.0_dp
    1355             :             END DO
    1356     3172680 :             point(1) = r(ir)*a(1, ia)
    1357     3172680 :             point(2) = r(ir)*a(2, ia)
    1358     3172680 :             point(3) = r(ir)*a(3, ia)
    1359    10785600 :             DO iatom = 1, natm_gauge
    1360     7612920 :                buf(iatom) = 1.0_dp
    1361    30451680 :                pra = point - ratom(:, iatom)
    1362     7612920 :                pa = SQRT(pra(1)**2 + pra(2)**2 + pra(3)**2)
    1363    31404240 :                DO jatom = 1, natm_gauge
    1364    20618640 :                   IF (iatom .EQ. jatom) CYCLE
    1365    52022880 :                   prb = point - ratom(:, jatom)
    1366    13005720 :                   pb = SQRT(prb(1)**2 + prb(2)**2 + prb(3)**2)
    1367    13005720 :                   ab = SQRT((pra(1) - prb(1))**2 + (pra(2) - prb(2))**2 + (pra(3) - prb(3))**2)
    1368             :                   !
    1369    13005720 :                   tmp = (pa - pb)/ab
    1370    13005720 :                   tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
    1371    13005720 :                   tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
    1372    13005720 :                   tmp = 0.5_dp*(3.0_dp - tmp**2)*tmp
    1373    28231560 :                   buf(iatom) = buf(iatom)*0.5_dp*(1.0_dp - tmp)
    1374             :                END DO
    1375             :             END DO
    1376    12755340 :             DO ixyz = 1, 3
    1377     9518040 :                res = 0.0_dp
    1378    32356800 :                DO iatom = 1, natm_gauge
    1379    32356800 :                   res = res + ratom(ixyz, iatom)*buf(iatom)
    1380             :                END DO
    1381    32356800 :                res = res/SUM(buf(1:natm_gauge))
    1382             :                !
    1383    12690720 :                g(ixyz, ir, ia) = res
    1384             :             END DO
    1385             :          END DO
    1386             :          END DO
    1387        1323 :          DEALLOCATE (buf)
    1388        1323 :       END SUBROUTINE get_gauge
    1389             :    END SUBROUTINE calculate_jrho_atom_ang
    1390             : 
    1391             : ! **************************************************************************************************
    1392             : !> \brief ...
    1393             : !> \param current_env ...
    1394             : !> \param qs_env ...
    1395             : !> \param iB ...
    1396             : !> \param idir ...
    1397             : ! **************************************************************************************************
    1398         864 :    SUBROUTINE calculate_jrho_atom(current_env, qs_env, iB, idir)
    1399             :       TYPE(current_env_type)                             :: current_env
    1400             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1401             :       INTEGER, INTENT(IN)                                :: iB, idir
    1402             : 
    1403             :       INTEGER                                            :: iat, iatom, ikind, ispin, jatom, &
    1404             :                                                             natm_gauge, natm_tot, natom, nkind, &
    1405             :                                                             nspins
    1406             :       INTEGER, DIMENSION(2)                              :: bo
    1407         864 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
    1408             :       LOGICAL                                            :: do_igaim, gapw, paw_atom
    1409             :       REAL(dp)                                           :: hard_radius, r(3)
    1410             :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: ratom
    1411         864 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1412             :       TYPE(cell_type), POINTER                           :: cell
    1413             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1414             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1415             :       TYPE(grid_atom_type), POINTER                      :: grid_atom
    1416             :       TYPE(harmonics_atom_type), POINTER                 :: harmonics
    1417         864 :       TYPE(jrho_atom_type), DIMENSION(:), POINTER        :: jrho1_atom_set
    1418             :       TYPE(jrho_atom_type), POINTER                      :: jrho1_atom
    1419         864 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1420         864 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1421             : 
    1422         864 :       NULLIFY (para_env, dft_control)
    1423         864 :       NULLIFY (jrho1_atom_set, grid_atom, harmonics)
    1424         864 :       NULLIFY (atomic_kind_set, qs_kind_set, atom_list)
    1425             : 
    1426             :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
    1427             :                       atomic_kind_set=atomic_kind_set, &
    1428             :                       qs_kind_set=qs_kind_set, &
    1429             :                       particle_set=particle_set, &
    1430             :                       cell=cell, &
    1431         864 :                       para_env=para_env)
    1432             : 
    1433             :       CALL get_current_env(current_env=current_env, &
    1434         864 :                            jrho1_atom_set=jrho1_atom_set)
    1435             : 
    1436         864 :       do_igaim = .FALSE.
    1437         864 :       IF (current_env%gauge .EQ. current_gauge_atom) do_igaim = .TRUE.
    1438             : 
    1439         864 :       gapw = dft_control%qs_control%gapw
    1440         864 :       nkind = SIZE(qs_kind_set, 1)
    1441         864 :       nspins = dft_control%nspins
    1442             : 
    1443         864 :       natm_tot = SIZE(particle_set)
    1444        2592 :       ALLOCATE (ratom(3, natm_tot))
    1445             : 
    1446         864 :       IF (gapw) THEN
    1447        2466 :          DO ikind = 1, nkind
    1448        1602 :             NULLIFY (atom_list, grid_atom, harmonics)
    1449        1602 :             CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
    1450             :             CALL get_qs_kind(qs_kind_set(ikind), &
    1451             :                              grid_atom=grid_atom, &
    1452             :                              harmonics=harmonics, &
    1453             :                              hard_radius=hard_radius, &
    1454        1602 :                              paw_atom=paw_atom)
    1455             : 
    1456        1602 :             IF (.NOT. paw_atom) CYCLE
    1457             : 
    1458             :             ! Distribute the atoms of this kind
    1459             : 
    1460        1242 :             bo = get_limit(natom, para_env%num_pe, para_env%mepos)
    1461             : 
    1462        4554 :             DO iat = bo(1), bo(2)
    1463         846 :                iatom = atom_list(iat)
    1464             :                NULLIFY (jrho1_atom)
    1465         846 :                jrho1_atom => jrho1_atom_set(iatom)
    1466             : 
    1467         846 :                natm_gauge = 0
    1468        3690 :                DO jatom = 1, natm_tot
    1469       11376 :                   r(:) = pbc(particle_set(jatom)%r(:) - particle_set(iatom)%r(:), cell)
    1470             :                   ! SQRT(SUM(r(:)**2)) .LE. 2.0_dp*hard_radius
    1471       12222 :                   IF (SUM(r(:)**2) .LE. (4.0_dp*hard_radius*hard_radius)) THEN
    1472        2160 :                      natm_gauge = natm_gauge + 1
    1473        8640 :                      ratom(:, natm_gauge) = r(:)
    1474             :                   END IF
    1475             :                END DO
    1476             : 
    1477        3771 :                DO ispin = 1, nspins
    1478     3237237 :                   jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef = 0.0_dp
    1479     3237237 :                   jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef = 0.0_dp
    1480             :                   CALL calculate_jrho_atom_ang(jrho1_atom, &
    1481             :                                                jrho1_atom%jrho_vec_rad_h(idir, ispin)%r_coef, &
    1482             :                                                jrho1_atom%jrho_vec_rad_s(idir, ispin)%r_coef, &
    1483             :                                                grid_atom, harmonics, &
    1484             :                                                do_igaim, &
    1485        2169 :                                                ratom, natm_gauge, iB, idir, ispin)
    1486             :                END DO !ispin
    1487             :             END DO !iat
    1488             :          END DO !ikind
    1489             :       END IF
    1490             : 
    1491         864 :       DEALLOCATE (ratom)
    1492             : 
    1493         864 :    END SUBROUTINE calculate_jrho_atom
    1494             : 
    1495             : END MODULE qs_linres_atom_current

Generated by: LCOV version 1.15