LCOV - code coverage report
Current view: top level - src/xc - xc.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 70.3 % 1278 898
Test Date: 2025-07-25 12:55:17 Functions: 93.3 % 30 28

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Exchange and Correlation functional calculations
      10              : !> \par History
      11              : !>      (13-Feb-2001) JGH, based on earlier version of apsi
      12              : !>      02.2003 Many many changes [fawzi]
      13              : !>      03.2004 new xc interface [fawzi]
      14              : !>      04.2004 kinetic functionals [fawzi]
      15              : !> \author fawzi
      16              : ! **************************************************************************************************
      17              : MODULE xc
      18              :    #:include 'xc.fypp'
      19              :    USE cp_array_utils, ONLY: cp_3d_r_cp_type
      20              :    USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next, &
      21              :                                       cp_sll_xc_deriv_type
      22              :    USE cp_log_handling, ONLY: cp_get_default_logger, &
      23              :                               cp_logger_get_default_unit_nr, &
      24              :                               cp_logger_type, &
      25              :                               cp_to_string
      26              :    USE input_section_types, ONLY: section_get_ival, &
      27              :                                   section_get_lval, &
      28              :                                   section_get_rval, &
      29              :                                   section_vals_get_subs_vals, &
      30              :                                   section_vals_type, &
      31              :                                   section_vals_val_get
      32              :    USE kahan_sum, ONLY: accurate_dot_product, &
      33              :                         accurate_sum
      34              :    USE kinds, ONLY: default_path_length, &
      35              :                     dp
      36              :    USE pw_grid_types, ONLY: PW_MODE_DISTRIBUTED, &
      37              :                             pw_grid_type
      38              :    USE pw_methods, ONLY: pw_axpy, &
      39              :                          pw_copy, &
      40              :                          pw_derive, &
      41              :                          pw_scale, &
      42              :                          pw_transfer, &
      43              :                          pw_zero, pw_integrate_function
      44              :    USE pw_pool_types, ONLY: &
      45              :       pw_pool_type
      46              :    USE pw_types, ONLY: &
      47              :       pw_c1d_gs_type, pw_r3d_rs_type
      48              :    USE xc_derivative_desc, ONLY: &
      49              :       deriv_rho, deriv_rhoa, deriv_rhob, &
      50              :       deriv_norm_drhoa, deriv_norm_drhob, deriv_norm_drho, deriv_tau_a, deriv_tau_b, deriv_tau, &
      51              :       deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, id_to_desc
      52              :    USE xc_derivative_set_types, ONLY: xc_derivative_set_type, &
      53              :                                       xc_dset_create, &
      54              :                                       xc_dset_get_derivative, &
      55              :                                       xc_dset_release, &
      56              :                                       xc_dset_zero_all, xc_dset_recover_pw
      57              :    USE xc_derivative_types, ONLY: xc_derivative_get, &
      58              :                                   xc_derivative_type
      59              :    USE xc_derivatives, ONLY: xc_functionals_eval, &
      60              :                              xc_functionals_get_needs
      61              :    USE xc_input_constants, ONLY: &
      62              :       xc_debug_new_routine, xc_default_f_routine, xc_test_lsd_f_routine
      63              :    USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
      64              :    USE xc_rho_set_types, ONLY: xc_rho_set_create, &
      65              :                                xc_rho_set_get, &
      66              :                                xc_rho_set_release, &
      67              :                                xc_rho_set_type, &
      68              :                                xc_rho_set_update, xc_rho_set_recover_pw
      69              :    USE xc_util, ONLY: xc_pw_smooth, xc_pw_laplace, xc_pw_divergence, xc_requires_tmp_g
      70              : #include "../base/base_uses.f90"
      71              : 
      72              :    IMPLICIT NONE
      73              :    PRIVATE
      74              :    PUBLIC :: xc_vxc_pw_create1, xc_vxc_pw_create, &
      75              :              xc_exc_calc, xc_calc_2nd_deriv_analytical, xc_calc_2nd_deriv_numerical, &
      76              :              xc_calc_2nd_deriv, xc_prep_2nd_deriv, divide_by_norm_drho, smooth_cutoff, &
      77              :              xc_uses_kinetic_energy_density, xc_uses_norm_drho
      78              :    PUBLIC :: calc_xc_density
      79              : 
      80              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      81              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc'
      82              : 
      83              : CONTAINS
      84              : 
      85              : ! **************************************************************************************************
      86              : !> \brief ...
      87              : !> \param xc_fun_section ...
      88              : !> \param lsd ...
      89              : !> \return ...
      90              : ! **************************************************************************************************
      91         7624 :    FUNCTION xc_uses_kinetic_energy_density(xc_fun_section, lsd) RESULT(res)
      92              :       TYPE(section_vals_type), POINTER, INTENT(IN) :: xc_fun_section
      93              :       LOGICAL, INTENT(IN) :: lsd
      94              :       LOGICAL :: res
      95              : 
      96              :       TYPE(xc_rho_cflags_type)                           :: needs
      97              : 
      98              :       needs = xc_functionals_get_needs(xc_fun_section, &
      99              :                                        lsd=lsd, &
     100        15248 :                                        calc_potential=.FALSE.)
     101         7624 :       res = (needs%tau_spin .OR. needs%tau)
     102              : 
     103         7624 :    END FUNCTION
     104              : 
     105              : ! **************************************************************************************************
     106              : !> \brief ...
     107              : !> \param xc_fun_section ...
     108              : !> \param lsd ...
     109              : !> \return ...
     110              : ! **************************************************************************************************
     111         7640 :    FUNCTION xc_uses_norm_drho(xc_fun_section, lsd) RESULT(res)
     112              :       TYPE(section_vals_type), POINTER, INTENT(IN) :: xc_fun_section
     113              :       LOGICAL, INTENT(IN) :: lsd
     114              :       LOGICAL :: res
     115              : 
     116              :       TYPE(xc_rho_cflags_type)                           :: needs
     117              : 
     118              :       needs = xc_functionals_get_needs(xc_fun_section, &
     119              :                                        lsd=lsd, &
     120        15280 :                                        calc_potential=.FALSE.)
     121         7640 :       res = (needs%norm_drho .OR. needs%norm_drho_spin)
     122              : 
     123         7640 :    END FUNCTION
     124              : 
     125              : ! **************************************************************************************************
     126              : !> \brief Exchange and Correlation functional calculations.
     127              : !>      depending on the selected functional_routine calls
     128              : !>      the correct routine
     129              : !> \param vxc_rho will contain the v_xc part that depend on rho
     130              : !>        (if one of the chosen xc functionals has it it is allocated and you
     131              : !>        are responsible for it)
     132              : !> \param vxc_tau will contain the kinetic tau part of v_xcthe functional
     133              : !>        (if one of the chosen xc functionals has it it is allocated and you
     134              : !>        are responsible for it)
     135              : !> \param rho_r the value of the density in the real space
     136              : !> \param rho_g value of the density in the g space (needs to be associated
     137              : !>        only for gradient corrections)
     138              : !> \param tau value of the kinetic density tau on the grid (can be null,
     139              : !>        used only with meta functionals)
     140              : !> \param exc the xc energy
     141              : !> \param xc_section parameters selecting the xc and the method used to
     142              : !>        calculate it
     143              : !> \param pw_pool the pool for the grids
     144              : !> \param compute_virial should the virial be computed... if so virial_xc must be present
     145              : !> \param virial_xc for calculating the GGA part of the stress
     146              : !> \param exc_r ...
     147              : !> \author fawzi
     148              : ! **************************************************************************************************
     149       218542 :    SUBROUTINE xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, &
     150              :                                 pw_pool, compute_virial, virial_xc, exc_r)
     151              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: vxc_rho, vxc_tau
     152              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau
     153              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
     154              :       REAL(KIND=dp), INTENT(out)                         :: exc
     155              :       TYPE(section_vals_type), POINTER                   :: xc_section
     156              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     157              :       LOGICAL, INTENT(IN)                                :: compute_virial
     158              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: virial_xc
     159              :       TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL           :: exc_r
     160              : 
     161              :       INTEGER                                            :: f_routine
     162              : 
     163       109271 :       CPASSERT(ASSOCIATED(rho_r))
     164       109271 :       CPASSERT(ASSOCIATED(xc_section))
     165       109271 :       CPASSERT(.NOT. ASSOCIATED(vxc_rho))
     166       109271 :       CPASSERT(.NOT. ASSOCIATED(vxc_tau))
     167              : 
     168       109271 :       virial_xc = 0.0_dp
     169              : 
     170              :       CALL section_vals_val_get(xc_section, "FUNCTIONAL_ROUTINE", &
     171       109271 :                                 i_val=f_routine)
     172       109271 :       SELECT CASE (f_routine)
     173              :       CASE (xc_default_f_routine)
     174              :          CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, tau=tau, exc_r=exc_r, &
     175              :                                rho_r=rho_r, rho_g=rho_g, exc=exc, xc_section=xc_section, &
     176       109271 :                                pw_pool=pw_pool, compute_virial=compute_virial, virial_xc=virial_xc)
     177              :       CASE (xc_debug_new_routine)
     178            0 :          CPASSERT(.NOT. PRESENT(exc_r))
     179              :          CALL xc_vxc_pw_create_debug(vxc_rho=vxc_rho, vxc_tau=vxc_tau, tau=tau, &
     180              :                                      rho_r=rho_r, rho_g=rho_g, exc=exc, xc_section=xc_section, &
     181            0 :                                      pw_pool=pw_pool)
     182              :       CASE (xc_test_lsd_f_routine)
     183            0 :          CPASSERT(.NOT. PRESENT(exc_r))
     184              :          CALL xc_vxc_pw_create_test_lsd(vxc_rho=vxc_rho, vxc_tau=vxc_tau, &
     185              :                                         tau=tau, rho_r=rho_r, rho_g=rho_g, exc=exc, &
     186       109271 :                                         xc_section=xc_section, pw_pool=pw_pool)
     187              :       CASE default
     188              :       END SELECT
     189              : 
     190       109271 :    END SUBROUTINE xc_vxc_pw_create1
     191              : 
     192              : ! **************************************************************************************************
     193              : !> \brief calculates vxc using lsd with rhoa=rhob=0.5*rho and compares
     194              : !>      with the lda result
     195              : !> \param vxc_rho will contain the v_xc part that depend on rho
     196              : !>        (if one of the chosen xc functionals has it it is allocated and you
     197              : !>        are responsible for it)
     198              : !> \param vxc_tau will contain the kinetic tau part of v_xc
     199              : !>        (if one of the chosen xc functionals has it it is allocated and you
     200              : !>        are responsible for it)
     201              : !> \param rho_r the value of the density in the real space
     202              : !> \param rho_g value of the density in the g space (needs to be associated
     203              : !>        only for gradient corrections)
     204              : !> \param tau value of the kinetic density tau on the grid (can be null,
     205              : !>        used only with meta functionals)
     206              : !> \param exc the xc energy
     207              : !> \param xc_section which functional to calculate, and how
     208              : !> \param pw_pool the pool for the grids
     209              : !> \author Fawzi Mohamed
     210              : !> \note
     211              : !>      for debugging only: leaks, and non parallel
     212              : ! **************************************************************************************************
     213            0 :    SUBROUTINE xc_vxc_pw_create_test_lsd(vxc_rho, vxc_tau, rho_r, rho_g, tau, &
     214              :                                         exc, xc_section, pw_pool)
     215              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: vxc_rho, vxc_tau
     216              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               :: rho_r, tau
     217              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
     218              :       REAL(KIND=dp), INTENT(out)                         :: exc
     219              :       TYPE(section_vals_type), POINTER                   :: xc_section
     220              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     221              : 
     222              :       LOGICAL, PARAMETER                                 :: compute_virial = .FALSE.
     223              : 
     224              :       CHARACTER(len=default_path_length)                 :: filename
     225            0 :       INTEGER, DIMENSION(:), POINTER         :: split_desc
     226              :       INTEGER                                            :: i, ii, ispin, j, k, order
     227              :       INTEGER, DIMENSION(2, 3)                           :: bo
     228              :       REAL(kind=dp)                                      :: diff, exc2, maxdiff, tmp
     229              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_xc
     230            0 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: pot, pot2, pot3
     231              :       TYPE(cp_sll_xc_deriv_type), POINTER                :: deriv_iter
     232            0 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: vxc_rho2, vxc_tau2
     233            0 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho2_g
     234            0 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               ::  rho2_r, tau2
     235              :       TYPE(xc_derivative_set_type)                       :: dSet1, dSet2
     236              :       TYPE(xc_derivative_type), POINTER                  :: deriv, deriv2, deriv3
     237              :       TYPE(xc_rho_set_type)                              :: rho_set1, rho_set2
     238              : 
     239            0 :       NULLIFY (vxc_rho2, vxc_tau2, tau2, pot, pot3, pot3, &
     240            0 :                deriv, deriv2, deriv3, rho2_g)
     241              : 
     242            0 :       bo = rho_r(1)%pw_grid%bounds_local
     243              : 
     244            0 :       ALLOCATE (rho2_r(2))
     245            0 :       DO ispin = 1, 2
     246            0 :          CALL pw_pool%create_pw(rho2_r(ispin))
     247              :       END DO
     248            0 :       DO k = bo(1, 3), bo(2, 3)
     249            0 :          DO j = bo(1, 2), bo(2, 2)
     250            0 :             DO i = bo(1, 1), bo(2, 1)
     251            0 :                tmp = rho_r(1)%array(i, j, k)*0.5
     252            0 :                rho2_r(1)%array(i, j, k) = tmp
     253            0 :                rho2_r(2)%array(i, j, k) = tmp
     254              :             END DO
     255              :          END DO
     256              :       END DO
     257              : 
     258            0 :       IF (ASSOCIATED(tau)) THEN
     259            0 :          ALLOCATE (tau2(2))
     260            0 :          DO ispin = 1, 2
     261            0 :             CALL pw_pool%create_pw(tau2(ispin))
     262              :          END DO
     263              : 
     264            0 :          DO k = bo(1, 3), bo(2, 3)
     265            0 :             DO j = bo(1, 2), bo(2, 2)
     266            0 :                DO i = bo(1, 1), bo(2, 1)
     267            0 :                   tmp = tau(1)%array(i, j, k)*0.5
     268            0 :                   tau2(1)%array(i, j, k) = tmp
     269            0 :                   tau2(2)%array(i, j, k) = tmp
     270              :                END DO
     271              :             END DO
     272              :          END DO
     273              :       END IF
     274              : 
     275            0 :       PRINT *, "about to calculate xc (lda)"
     276              :       CALL xc_rho_set_and_dset_create(rho_r=rho_r, rho_g=rho_g, &
     277              :                                       tau=tau, xc_section=xc_section, &
     278              :                                       pw_pool=pw_pool, rho_set=rho_set1, &
     279              :                                       deriv_set=dSet1, deriv_order=1, &
     280            0 :                                       calc_potential=.FALSE.)
     281              :       CALL xc_vxc_pw_create(rho_r=rho_r, rho_g=rho_g, tau=tau, &
     282              :                             vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
     283              :                             pw_pool=pw_pool, &
     284            0 :                             compute_virial=compute_virial, virial_xc=virial_xc)
     285            0 :       PRINT *, "did calculate xc (lda)"
     286            0 :       PRINT *, "about to calculate xc (lsd)"
     287              :       CALL xc_rho_set_and_dset_create(rho_set=rho_set2, deriv_set=dSet2, &
     288              :                                       rho_r=rho2_r, rho_g=rho2_g, tau=tau2, xc_section=xc_section, &
     289              :                                       pw_pool=pw_pool, deriv_order=1, &
     290            0 :                                       calc_potential=.FALSE.)
     291              :       CALL xc_vxc_pw_create(rho_r=rho2_r, rho_g=rho2_g, tau=tau2, &
     292              :                             vxc_rho=vxc_rho2, vxc_tau=vxc_tau2, exc=exc2, xc_section=xc_section, &
     293            0 :                             pw_pool=pw_pool, compute_virial=compute_virial, virial_xc=virial_xc)
     294            0 :       PRINT *, "did calculate xc (new)"
     295            0 :       PRINT *, "at (0,0,0) rho_r=", rho_r(1)%array(0, 0, 0), &
     296            0 :          "rho2_r(1)=", rho2_r(1)%array(0, 0, 0), &
     297            0 :          "rho2_r(2)=", rho2_r(2)%array(0, 0, 0), &
     298            0 :          "rho_r_sm=", rho_set1%rho(0, 0, 0), "rhoa2_r_sm=", rho_set2%rhoa(0, 0, 0), &
     299            0 :          "rhob2_r_sm=", rho_set2%rhob(0, 0, 0)
     300              :       OPEN (unit=120, file="rho.bindata", status="unknown", access='sequential', &
     301            0 :             form="unformatted", action="write")
     302            0 :       pot => rho_set1%rho
     303            0 :       WRITE (unit=120) pot(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(2, 3))
     304            0 :       CLOSE (unit=120)
     305              :       OPEN (unit=120, file="rhoa.bindata", status="unknown", access='sequential', &
     306            0 :             form="unformatted", action="write")
     307            0 :       pot => rho_set2%rhoa
     308            0 :       WRITE (unit=120) pot(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(2, 3))
     309            0 :       CLOSE (unit=120)
     310              :       OPEN (unit=120, file="rhob.bindata", status="unknown", access='sequential', &
     311            0 :             form="unformatted", action="write")
     312            0 :       pot => rho_set2%rhob
     313            0 :       WRITE (unit=120) pot(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(2, 3))
     314            0 :       CLOSE (unit=120)
     315              :       OPEN (unit=120, file="ndrho.bindata", status="unknown", access='sequential', &
     316            0 :             form="unformatted", action="write")
     317            0 :       pot => rho_set1%norm_drho
     318            0 :       WRITE (unit=120) pot(:, :, bo(2, 3))
     319            0 :       CLOSE (unit=120)
     320              :       OPEN (unit=120, file="ndrhoa.bindata", status="unknown", access='sequential', &
     321            0 :             form="unformatted", action="write")
     322            0 :       pot => rho_set2%norm_drhoa
     323            0 :       WRITE (unit=120) pot(:, :, bo(2, 3))
     324            0 :       CLOSE (unit=120)
     325              :       OPEN (unit=120, file="ndrhob.bindata", status="unknown", access='sequential', &
     326            0 :             form="unformatted", action="write")
     327            0 :       pot => rho_set2%norm_drhob
     328            0 :       WRITE (unit=120) pot(:, :, bo(2, 3))
     329            0 :       CLOSE (unit=120)
     330            0 :       IF (rho_set1%has%tau) THEN
     331              :          OPEN (unit=120, file="tau.bindata", status="unknown", access='sequential', &
     332            0 :                form="unformatted", action="write")
     333            0 :          pot => rho_set1%tau
     334            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     335            0 :          CLOSE (unit=120)
     336              :       END IF
     337            0 :       IF (rho_set2%has%tau_spin) THEN
     338              :          OPEN (unit=120, file="tau_a.bindata", status="unknown", access='sequential', &
     339            0 :                form="unformatted", action="write")
     340            0 :          pot => rho_set2%tau_a
     341            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     342            0 :          CLOSE (unit=120)
     343              :          OPEN (unit=120, file="tau_v.bindata", status="unknown", access='sequential', &
     344            0 :                form="unformatted", action="write")
     345            0 :          pot => rho_set2%tau_b
     346            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     347            0 :          CLOSE (unit=120)
     348              :       END IF
     349              :       OPEN (unit=120, file="vxc.bindata", status="unknown", access='sequential', &
     350            0 :             form="unformatted", action="write")
     351            0 :       pot => vxc_rho(1)%array
     352            0 :       WRITE (unit=120) pot(:, :, bo(2, 3))
     353            0 :       CLOSE (unit=120)
     354              :       OPEN (unit=120, file="vxc2.bindata", status="unknown", access='sequential', &
     355            0 :             form="unformatted", action="write")
     356            0 :       pot => vxc_rho2(1)%array
     357            0 :       WRITE (unit=120) pot(:, :, bo(2, 3))
     358            0 :       CLOSE (unit=120)
     359            0 :       IF (ASSOCIATED(vxc_tau)) THEN
     360              :          OPEN (unit=120, file="vxc_tau.bindata", status="unknown", access='sequential', &
     361            0 :                form="unformatted", action="write")
     362            0 :          pot => vxc_tau(1)%array
     363            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     364            0 :          CLOSE (unit=120)
     365              :       END IF
     366            0 :       IF (ASSOCIATED(vxc_tau2)) THEN
     367              :          OPEN (unit=120, file="vxc_tau2_a.bindata", status="unknown", access='sequential', &
     368            0 :                form="unformatted", action="write")
     369            0 :          pot => vxc_tau2(1)%array
     370            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     371            0 :          CLOSE (unit=120)
     372              :          OPEN (unit=120, file="vxc_tau2_b.bindata", status="unknown", access='sequential', &
     373            0 :                form="unformatted", action="write")
     374            0 :          pot => vxc_tau2(2)%array
     375            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     376            0 :          CLOSE (unit=120)
     377              :       END IF
     378              : 
     379            0 :       PRINT *, "calc diff on vxc"
     380            0 :       maxDiff = 0.0_dp
     381            0 :       DO ispin = 1, 1
     382            0 :          ii = 0
     383            0 :          DO k = bo(1, 3), bo(2, 3)
     384            0 :             DO j = bo(1, 2), bo(2, 2)
     385            0 :                DO i = bo(1, 1), bo(2, 1)
     386            0 :                   ii = ii + 1
     387              :                   diff = ABS(vxc_rho(ispin)%array(i, j, k) - &
     388            0 :                              vxc_rho2(ispin)%array(i, j, k))
     389            0 :                   IF (ii == 1) THEN
     390            0 :                      PRINT *, "vxc", ispin, "=", vxc_rho(ispin)%array(i, j, k), &
     391            0 :                         "vs", vxc_rho2(ispin)%array(i, j, k), "diff=", diff
     392              :                   END IF
     393            0 :                   IF (maxDiff < diff) THEN
     394            0 :                      maxDiff = diff
     395            0 :                      PRINT *, "diff=", diff, " at ", i, ",", j, ",", k, &
     396            0 :                         " spin=", ispin, "rho=", rho_set1%rho(i, j, k), &
     397            0 :                         " ndrho=", rho_set1%norm_drho(i, j, k)
     398              :                   END IF
     399              :                END DO
     400              :             END DO
     401              :          END DO
     402              :       END DO
     403            0 :       PRINT *, "diff exc=", ABS(exc - exc2), "diff vxc=", maxdiff
     404              :       ! CPASSERT(maxdiff<5.e-11)
     405              :       ! CPASSERT(ABS(exc-exc2)<1.e-14)
     406              : 
     407            0 :       IF (ASSOCIATED(vxc_tau)) THEN
     408            0 :          PRINT *, "calc diff on vxc_tau"
     409            0 :          maxDiff = 0.0_dp
     410            0 :          DO ispin = 1, 1
     411            0 :             ii = 0
     412            0 :             DO k = bo(1, 3), bo(2, 3)
     413            0 :                DO j = bo(1, 2), bo(2, 2)
     414            0 :                   DO i = bo(1, 1), bo(2, 1)
     415            0 :                      ii = ii + 1
     416              :                      diff = ABS(vxc_tau(ispin)%array(i, j, k) - &
     417            0 :                                 vxc_tau2(ispin)%array(i, j, k))
     418            0 :                      IF (ii == 1) THEN
     419            0 :                         PRINT *, "vxc_tau", ispin, "=", vxc_tau(ispin)%array(i, j, k), &
     420            0 :                            "vs", vxc_tau2(ispin)%array(i, j, k), "diff=", diff
     421              :                      END IF
     422            0 :                      IF (maxDiff < diff) THEN
     423            0 :                         maxDiff = diff
     424            0 :                         PRINT *, "diff=", diff, " at ", i, ",", j, ",", k, &
     425            0 :                            " spin=", ispin, "rho=", rho_set1%rho(i, j, k), &
     426            0 :                            " ndrho=", rho_set1%norm_drho(i, j, k)
     427              :                      END IF
     428              :                   END DO
     429              :                END DO
     430              :             END DO
     431              :          END DO
     432            0 :          PRINT *, "diff exc=", ABS(exc - exc2), "diff vxc_tau=", maxdiff
     433              :       END IF
     434            0 :       deriv_iter => dSet1%derivs
     435            0 :       DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
     436              :          CALL xc_derivative_get(deriv, order=order, &
     437            0 :                                 split_desc=split_desc, deriv_data=pot)
     438            0 :          SELECT CASE (order)
     439              :          CASE (0)
     440            0 :             filename = "e_0.bindata"
     441            0 :             deriv2 => xc_dset_get_derivative(dSet2, [INTEGER::])
     442              :          CASE (1)
     443            0 :             filename = "e_"//TRIM(id_to_desc(split_desc(1)))//".bindata"
     444            0 :             IF (split_desc(1) == deriv_rho) THEN
     445            0 :                deriv2 => xc_dset_get_derivative(dSet2, [deriv_rhoa])
     446            0 :             ELSEIF (split_desc(1) == deriv_tau) THEN
     447            0 :                deriv2 => xc_dset_get_derivative(dSet2, [deriv_tau_a])
     448            0 :             ELSEIF (split_desc(1) == deriv_norm_drho) THEN
     449            0 :                deriv2 => xc_dset_get_derivative(dSet2, [deriv_norm_drhoa])
     450            0 :                deriv3 => xc_dset_get_derivative(dSet2, [deriv_norm_drho])
     451            0 :                IF (ASSOCIATED(deriv3)) THEN
     452            0 :                   IF (ASSOCIATED(deriv2)) THEN
     453              :                      CALL xc_derivative_get(deriv2, &
     454            0 :                                             deriv_data=pot2)
     455              :                      CALL xc_derivative_get(deriv3, &
     456            0 :                                             deriv_data=pot3)
     457            0 :                      pot2 = pot2 + pot3
     458              :                   ELSE
     459              :                      deriv2 => deriv3
     460              :                   END IF
     461            0 :                   NULLIFY (deriv3, pot2, pot3)
     462              :                END IF
     463              :             ELSE
     464            0 :                CPABORT("Unknown derivative variable")
     465              :             END IF
     466              :          CASE default
     467            0 :             CPABORT("Unsupported derivative order")
     468              :          END SELECT
     469              :          CALL xc_derivative_get(deriv2, &
     470            0 :                                 deriv_data=pot2)
     471            0 :          PRINT *, "checking ", filename
     472            0 :          maxDiff = 0.0_dp
     473            0 :          DO k = bo(1, 3), bo(2, 3)
     474            0 :             DO j = bo(1, 2), bo(2, 2)
     475            0 :                DO i = bo(1, 1), bo(2, 1)
     476            0 :                   diff = ABS(pot(i, j, k) - pot2(i, j, k))
     477            0 :                   IF (maxDiff < diff) THEN
     478            0 :                      maxDiff = diff
     479            0 :                      PRINT *, "ediff(", i, j, k, ")=", maxDiff, &
     480            0 :                         "rho=", rho_set1%rho(i, j, k), &
     481            0 :                         "ndrho=", rho_set1%norm_drho(i, j, k)
     482              :                   END IF
     483              :                END DO
     484              :             END DO
     485              :          END DO
     486            0 :          PRINT *, "maxdiff ", filename, "=", maxDiff
     487              :          OPEN (unit=120, file=TRIM(filename), status="unknown", &
     488              :                access='sequential', &
     489            0 :                form="unformatted")
     490            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     491            0 :          CLOSE (unit=120)
     492              :       END DO
     493            0 :       deriv_iter => dSet2%derivs
     494            0 :       DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
     495              :          CALL xc_derivative_get(deriv, order=order, &
     496            0 :                                 split_desc=split_desc, deriv_data=pot)
     497            0 :          SELECT CASE (order)
     498              :          CASE (0)
     499            0 :             filename = "e_0-2.bindata"
     500              :          CASE (1)
     501            0 :             filename = "e_"//TRIM(id_to_desc(split_desc(1)))//"-2.bindata"
     502              :          CASE default
     503            0 :             CPABORT("Unsupported derivative order")
     504              :          END SELECT
     505              :          OPEN (unit=120, file=TRIM(filename), status="unknown", &
     506              :                access='sequential', &
     507            0 :                form="unformatted")
     508            0 :          WRITE (unit=120) pot(:, :, bo(2, 3))
     509            0 :          CLOSE (unit=120)
     510              :       END DO
     511            0 :       CALL xc_rho_set_release(rho_set1)
     512            0 :       CALL xc_rho_set_release(rho_set2)
     513            0 :       CALL xc_dset_release(dSet2)
     514            0 :       CALL xc_dset_release(dSet1)
     515            0 :       DO ispin = 1, 2
     516            0 :          CALL pw_pool%give_back_pw(rho2_r(ispin))
     517            0 :          CALL pw_pool%give_back_pw(vxc_rho2(ispin))
     518            0 :          IF (ASSOCIATED(vxc_tau2)) THEN
     519            0 :             CALL pw_pool%give_back_pw(vxc_tau2(ispin))
     520              :          END IF
     521              :       END DO
     522            0 :       DEALLOCATE (vxc_rho2, rho2_r, rho2_g)
     523            0 :       IF (ASSOCIATED(vxc_tau2)) THEN
     524            0 :          DEALLOCATE (vxc_tau2)
     525              :       END IF
     526              : 
     527            0 :    END SUBROUTINE xc_vxc_pw_create_test_lsd
     528              : 
     529              : ! **************************************************************************************************
     530              : !> \brief calculates vxc outputting the yz plane of rho, and of the various components
     531              : !>      of the the derivatives and of vxc
     532              : !> \param vxc_rho will contain the v_xc part that depend on rho
     533              : !>        (if one of the chosen xc functionals has it it is allocated and you
     534              : !>        are responsible for it)
     535              : !> \param vxc_tau will contain the kinetic tau part of v_xc
     536              : !>        (if one of the chosen xc functionals has it it is allocated and you
     537              : !>        are responsible for it)
     538              : !> \param rho_r the value of the density in the real space
     539              : !> \param rho_g value of the density in the g space (needs to be associated
     540              : !>        only for gradient corrections)
     541              : !> \param tau value of the kinetic density tau on the grid (can be null,
     542              : !>        used only with meta functionals)
     543              : !> \param exc the xc energy
     544              : !> \param xc_section which functional should be used, and how to do it
     545              : !> \param pw_pool the pool for the grids
     546              : !> \author Fawzi Mohamed
     547              : !> \note
     548              : !>      for debugging only.
     549              : ! **************************************************************************************************
     550            0 :    SUBROUTINE xc_vxc_pw_create_debug(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, &
     551              :                                      xc_section, pw_pool)
     552              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: vxc_rho, vxc_tau
     553              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               :: rho_r, tau
     554              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
     555              :       REAL(KIND=dp), INTENT(out)                         :: exc
     556              :       TYPE(section_vals_type), POINTER                   :: xc_section
     557              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     558              : 
     559              :       LOGICAL, PARAMETER                                 :: compute_virial = .FALSE.
     560              : 
     561              :       CHARACTER(len=default_path_length)                 :: filename
     562            0 :       INTEGER, DIMENSION(:), POINTER         :: split_desc
     563              :       INTEGER                                            :: i, ispin, j, k, order
     564              :       INTEGER, DIMENSION(2, 3)                           :: bo
     565              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_xc
     566            0 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: pot
     567              :       TYPE(cp_logger_type), POINTER                      :: logger
     568              :       TYPE(cp_sll_xc_deriv_type), POINTER                :: deriv_iter
     569              :       TYPE(xc_derivative_set_type)                       :: dSet1
     570              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     571              :       TYPE(xc_rho_set_type)                              :: rho_set1
     572              : 
     573            0 :       NULLIFY (pot, deriv)
     574            0 :       logger => cp_get_default_logger()
     575              : 
     576            0 :       bo = rho_r(1)%pw_grid%bounds_local
     577              : 
     578              :       CALL xc_rho_set_and_dset_create(rho_r=rho_r, rho_g=rho_g, &
     579              :                                       tau=tau, xc_section=xc_section, &
     580              :                                       pw_pool=pw_pool, rho_set=rho_set1, &
     581              :                                       deriv_set=dSet1, deriv_order=1, &
     582            0 :                                       calc_potential=.FALSE.)
     583              : 
     584              :       ! outputs 0,:,: plane
     585            0 :       IF (bo(1, 1) <= 0 .AND. 0 <= bo(2, 1)) THEN
     586            0 :          IF (rho_set1%has%rho_spin) THEN
     587              :             OPEN (unit=120, file="rhoa.bindata", status="unknown", access='sequential', &
     588            0 :                   form="unformatted", action="write")
     589            0 :             pot => rho_set1%rhoa
     590            0 :             WRITE (unit=120) pot(0, :, :)
     591            0 :             CLOSE (unit=120)
     592              :             OPEN (unit=120, file="rhob.bindata", status="unknown", access='sequential', &
     593            0 :                   form="unformatted", action="write")
     594            0 :             pot => rho_set1%rhob
     595            0 :             WRITE (unit=120) pot(0, :, :)
     596            0 :             CLOSE (unit=120)
     597              :          END IF
     598            0 :          IF (rho_set1%has%norm_drho) THEN
     599              :             OPEN (unit=120, file="ndrho.bindata", status="unknown", access='sequential', &
     600            0 :                   form="unformatted", action="write")
     601            0 :             pot => rho_set1%norm_drho
     602            0 :             WRITE (unit=120) pot(0, :, :)
     603            0 :             CLOSE (unit=120)
     604              :          END IF
     605            0 :          IF (rho_set1%has%norm_drho_spin) THEN
     606              :             OPEN (unit=120, file="ndrhoa.bindata", status="unknown", access='sequential', &
     607            0 :                   form="unformatted", action="write")
     608            0 :             pot => rho_set1%norm_drhoa
     609            0 :             WRITE (unit=120) pot(0, :, :)
     610            0 :             CLOSE (unit=120)
     611              :             OPEN (unit=120, file="ndrhob.bindata", status="unknown", access='sequential', &
     612            0 :                   form="unformatted", action="write")
     613            0 :             pot => rho_set1%norm_drhob
     614            0 :             WRITE (unit=120) pot(0, :, :)
     615            0 :             CLOSE (unit=120)
     616              :          END IF
     617            0 :          IF (rho_set1%has%rho) THEN
     618              :             OPEN (unit=120, file="rho.bindata", status="unknown", access='sequential', &
     619            0 :                   form="unformatted", action="write")
     620            0 :             pot => rho_set1%rho
     621            0 :             WRITE (unit=120) pot(0, :, :)
     622            0 :             CLOSE (unit=120)
     623              :          END IF
     624            0 :          IF (rho_set1%has%tau) THEN
     625              :             OPEN (unit=120, file="tau.bindata", status="unknown", access='sequential', &
     626            0 :                   form="unformatted", action="write")
     627            0 :             pot => rho_set1%tau
     628            0 :             WRITE (unit=120) pot(0, :, :)
     629            0 :             CLOSE (unit=120)
     630              :          END IF
     631            0 :          IF (rho_set1%has%tau_spin) THEN
     632              :             OPEN (unit=120, file="tau_a.bindata", status="unknown", access='sequential', &
     633            0 :                   form="unformatted", action="write")
     634            0 :             pot => rho_set1%tau_a
     635            0 :             WRITE (unit=120) pot(0, :, :)
     636            0 :             CLOSE (unit=120)
     637              :             OPEN (unit=120, file="tau_b.bindata", status="unknown", access='sequential', &
     638            0 :                   form="unformatted", action="write")
     639            0 :             pot => rho_set1%tau_b
     640            0 :             WRITE (unit=120) pot(0, :, :)
     641            0 :             CLOSE (unit=120)
     642              :          END IF
     643              : 
     644            0 :          deriv_iter => dSet1%derivs
     645            0 :          DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
     646              :             CALL xc_derivative_get(deriv, order=order, &
     647            0 :                                    split_desc=split_desc, deriv_data=pot)
     648            0 :             SELECT CASE (order)
     649              :             CASE (0)
     650            0 :                filename = "e_0.bindata"
     651              :             CASE (1)
     652            0 :                filename = "e_"//TRIM(id_to_desc(split_desc(1)))//".bindata"
     653              :             CASE default
     654            0 :                CPABORT("Unsupported derivative order")
     655              :             END SELECT
     656              :             OPEN (unit=120, file=TRIM(filename), status="unknown", &
     657              :                   access='sequential', &
     658            0 :                   form="unformatted")
     659            0 :             WRITE (unit=120) pot(0, :, :)
     660            0 :             CLOSE (unit=120)
     661              :          END DO
     662              :       END IF
     663              : 
     664              :       CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, &
     665              :                             rho_r=rho_r, rho_g=rho_g, tau=tau, &
     666              :                             exc=exc, xc_section=xc_section, &
     667              :                             pw_pool=pw_pool, &
     668            0 :                             compute_virial=compute_virial, virial_xc=virial_xc)
     669              : 
     670              :       ! outputs 0,:,: plane
     671            0 :       IF (bo(1, 1) <= 0 .AND. 0 <= bo(2, 1)) THEN
     672            0 :          IF (ASSOCIATED(vxc_rho)) THEN
     673            0 :             DO ispin = 1, SIZE(vxc_rho)
     674            0 :                WRITE (filename, "('vxc-',i1,'.bindata')") ispin
     675              :                OPEN (unit=120, file=filename, status="unknown", access='sequential', &
     676            0 :                      form="unformatted", action="write")
     677            0 :                pot => vxc_rho(ispin)%array
     678            0 :                WRITE (unit=120) pot(0, :, :)
     679            0 :                CLOSE (unit=120)
     680              : 
     681            0 :                pot => vxc_rho(ispin)%array
     682            0 :                DO k = bo(1, 3), bo(2, 3)
     683            0 :                   DO j = bo(1, 2), bo(2, 2)
     684            0 :                      DO i = bo(1, 1), bo(2, 1)
     685            0 :                         IF (ABS(pot(i, j, k)) > 10.0_dp) THEN
     686              :                            WRITE (cp_logger_get_default_unit_nr(logger), &
     687              :                                   "(' vxc_',i1,'(',i6,',',i6,',',i6,')=',e11.4)", &
     688            0 :                                   advance="no") ispin, i, j, k, pot(i, j, k)
     689            0 :                            IF (rho_set1%has%rho_spin) THEN
     690              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     691              :                                      "(' rho=(',e11.4,',',e11.4,')')", advance="no") &
     692            0 :                                  rho_set1%rhoa(i, j, k), rho_set1%rhob(i, j, k)
     693            0 :                            ELSE IF (rho_set1%has%rho) THEN
     694              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     695            0 :                                      "(' rho=',e11.4)", advance="no") rho_set1%rho(i, j, k)
     696              :                            END IF
     697            0 :                            IF (rho_set1%has%norm_drho_spin) THEN
     698              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     699              :                                      "(' ndrho=(',e11.4,',',e11.4,')')", advance="no") &
     700            0 :                                  rho_set1%norm_drhoa(i, j, k), &
     701            0 :                                  rho_set1%norm_drhob(i, j, k)
     702            0 :                            ELSE IF (rho_set1%has%norm_drho) THEN
     703              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     704            0 :                                      "(' ndrho=',e11.4)", advance="no") rho_set1%norm_drho(i, j, k)
     705              :                            END IF
     706            0 :                            IF (rho_set1%has%tau_spin) THEN
     707              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     708              :                                      "(' tau=(',e11.4,',',e11.4,')')", advance="no") &
     709            0 :                                  rho_set1%tau_a(i, j, k), &
     710            0 :                                  rho_set1%tau_b(i, j, k)
     711            0 :                            ELSE IF (rho_set1%has%tau) THEN
     712              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     713            0 :                                      "(' tau=',e11.4)", advance="no") rho_set1%tau(i, j, k)
     714              :                            END IF
     715              : 
     716            0 :                            deriv_iter => dSet1%derivs
     717            0 :                            DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
     718              :                               CALL xc_derivative_get(deriv, order=order, &
     719            0 :                                                      split_desc=split_desc, deriv_data=pot)
     720            0 :                               SELECT CASE (order)
     721              :                               CASE (0)
     722            0 :                                  filename = " e_0"
     723              :                               CASE (1)
     724            0 :                                  filename = " e_"//TRIM(id_to_desc(split_desc(1)))
     725              :                               CASE default
     726            0 :                                  CPABORT("Unsupported derivative order")
     727              :                               END SELECT
     728              :                               WRITE (unit=cp_logger_get_default_unit_nr(logger), &
     729              :                                      fmt="(a,'=',e11.4)", advance="no") &
     730            0 :                                  TRIM(filename), pot(i, j, k)
     731              :                            END DO
     732              : 
     733              :                            WRITE (cp_logger_get_default_unit_nr(logger), &
     734            0 :                                   "()")
     735              :                         END IF
     736              :                      END DO
     737              :                   END DO
     738              :                END DO
     739              :             END DO
     740              :          END IF
     741            0 :          IF (ASSOCIATED(vxc_tau)) THEN
     742            0 :             DO ispin = 1, SIZE(vxc_tau)
     743            0 :                WRITE (filename, "('vxc_tau_',i1,'.bindata')") ispin
     744              :                OPEN (unit=120, file=filename, status="unknown", access='sequential', &
     745            0 :                      form="unformatted", action="write")
     746            0 :                pot => vxc_tau(ispin)%array
     747            0 :                WRITE (unit=120) pot(0, :, :)
     748            0 :                CLOSE (unit=120)
     749              : 
     750            0 :                pot => vxc_tau(ispin)%array
     751            0 :                DO k = bo(1, 3), bo(2, 3)
     752            0 :                   DO j = bo(1, 2), bo(2, 2)
     753            0 :                      DO i = bo(1, 1), bo(2, 1)
     754            0 :                         IF (ABS(pot(i, j, k)) > 10.0_dp) THEN
     755              :                            WRITE (cp_logger_get_default_unit_nr(logger), &
     756              :                                   "('vxc_tau_',i1,'(',i6,',',i6,',',i6,')=',e11.4)", &
     757            0 :                                   advance="no") ispin, i, j, k, pot(i, j, k)
     758            0 :                            IF (rho_set1%has%rho_spin) THEN
     759              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     760              :                                      "(' rho=(',e11.4,',',e11.4,')')", advance="no") &
     761            0 :                                  rho_set1%rhoa(i, j, k), rho_set1%rhob(i, j, k)
     762            0 :                            ELSE IF (rho_set1%has%rho) THEN
     763              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     764            0 :                                      "(' rho=',e11.4)", advance="no") rho_set1%rho(i, j, k)
     765              :                            END IF
     766            0 :                            IF (rho_set1%has%norm_drho_spin) THEN
     767              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     768              :                                      "(' ndrho=(',e11.4,',',e11.4,')')", advance="no") &
     769            0 :                                  rho_set1%norm_drhoa(i, j, k), &
     770            0 :                                  rho_set1%norm_drhob(i, j, k)
     771            0 :                            ELSE IF (rho_set1%has%norm_drho) THEN
     772              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     773            0 :                                      "(' ndrho=',e11.4)", advance="no") rho_set1%norm_drho(i, j, k)
     774              :                            END IF
     775            0 :                            IF (rho_set1%has%tau_spin) THEN
     776              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     777              :                                      "(' tau=(',e11.4,',',e11.4,')')", advance="no") &
     778            0 :                                  rho_set1%tau_a(i, j, k), &
     779            0 :                                  rho_set1%tau_b(i, j, k)
     780            0 :                            ELSE IF (rho_set1%has%tau) THEN
     781              :                               WRITE (cp_logger_get_default_unit_nr(logger), &
     782            0 :                                      "(' tau=',e11.4)", advance="no") rho_set1%tau(i, j, k)
     783              :                            END IF
     784              : 
     785            0 :                            deriv_iter => dSet1%derivs
     786            0 :                            DO WHILE (cp_sll_xc_deriv_next(deriv_iter, el_att=deriv))
     787              :                               CALL xc_derivative_get(deriv, order=order, &
     788            0 :                                                      split_desc=split_desc, deriv_data=pot)
     789            0 :                               SELECT CASE (order)
     790              :                               CASE (0)
     791            0 :                                  filename = " e_0"
     792              :                               CASE (1)
     793            0 :                                  filename = " e_"//TRIM(id_to_desc(split_desc(1)))
     794              :                               CASE default
     795            0 :                                  CPABORT("Unsupported derivative order")
     796              :                               END SELECT
     797              :                               WRITE (unit=cp_logger_get_default_unit_nr(logger), &
     798              :                                      fmt="(a,'=',e11.4)", advance="no") &
     799            0 :                                  TRIM(filename), pot(i, j, k)
     800              :                            END DO
     801              : 
     802            0 :                            WRITE (cp_logger_get_default_unit_nr(logger), "()")
     803              :                         END IF
     804              :                      END DO
     805              :                   END DO
     806              :                END DO
     807              :             END DO
     808              :          END IF
     809              :       END IF
     810              : 
     811            0 :       CALL xc_dset_release(dSet1)
     812            0 :       CALL xc_rho_set_release(rho_set1)
     813            0 :    END SUBROUTINE xc_vxc_pw_create_debug
     814              : 
     815              : ! **************************************************************************************************
     816              : !> \brief creates a xc_rho_set and a derivative set containing the derivatives
     817              : !>      of the functionals with the given deriv_order.
     818              : !> \param rho_set will contain the rho set
     819              : !> \param deriv_set will contain the derivatives
     820              : !> \param deriv_order the order of the requested derivatives. If positive
     821              : !>        0:deriv_order are calculated, if negative only -deriv_order is
     822              : !>        guaranteed to be valid. Orders not requested might be present,
     823              : !>        but might contain garbage.
     824              : !> \param rho_r the value of the density in the real space
     825              : !> \param rho_g value of the density in the g space (can be null, used only
     826              : !>        without smoothing of rho or deriv)
     827              : !> \param tau value of the kinetic density tau on the grid (can be null,
     828              : !>        used only with meta functionals)
     829              : !> \param xc_section the section describing the functional to use
     830              : !> \param pw_pool the pool for the grids
     831              : !> \param calc_potential if the basic components of the arguments
     832              : !>        should be kept in rho set (a basic component is for example drho
     833              : !>        when with lda a functional needs norm_drho)
     834              : !> \author fawzi
     835              : !> \note
     836              : !>      if any of the functionals is gradient corrected the full gradient is
     837              : !>      added to the rho set
     838              : ! **************************************************************************************************
     839       267878 :    SUBROUTINE xc_rho_set_and_dset_create(rho_set, deriv_set, deriv_order, &
     840              :                                          rho_r, rho_g, tau, xc_section, pw_pool, &
     841              :                                          calc_potential)
     842              : 
     843              :       TYPE(xc_rho_set_type)                              :: rho_set
     844              :       TYPE(xc_derivative_set_type)                       :: deriv_set
     845              :       INTEGER, INTENT(in)                                :: deriv_order
     846              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               :: rho_r, tau
     847              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
     848              :       TYPE(section_vals_type), POINTER                   :: xc_section
     849              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     850              :       LOGICAL, INTENT(in)                                :: calc_potential
     851              : 
     852              :       CHARACTER(len=*), PARAMETER :: routineN = 'xc_rho_set_and_dset_create'
     853              : 
     854              :       INTEGER                                            :: handle, nspins
     855              :       LOGICAL                                            :: lsd
     856              :       TYPE(section_vals_type), POINTER                   :: xc_fun_sections
     857              : 
     858       133939 :       CALL timeset(routineN, handle)
     859              : 
     860       133939 :       CPASSERT(ASSOCIATED(pw_pool))
     861              : 
     862       133939 :       nspins = SIZE(rho_r)
     863       133939 :       lsd = (nspins /= 1)
     864              : 
     865       133939 :       xc_fun_sections => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     866              : 
     867       133939 :       CALL xc_dset_create(deriv_set, pw_pool)
     868              : 
     869              :       CALL xc_rho_set_create(rho_set, &
     870              :                              rho_r(1)%pw_grid%bounds_local, &
     871              :                              rho_cutoff=section_get_rval(xc_section, "density_cutoff"), &
     872              :                              drho_cutoff=section_get_rval(xc_section, "gradient_cutoff"), &
     873       133939 :                              tau_cutoff=section_get_rval(xc_section, "tau_cutoff"))
     874              : 
     875              :       CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, &
     876              :                              xc_functionals_get_needs(xc_fun_sections, lsd, calc_potential), &
     877              :                              section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
     878              :                              section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
     879       133939 :                              pw_pool)
     880              : 
     881              :       CALL xc_functionals_eval(xc_fun_sections, &
     882              :                                lsd=lsd, &
     883              :                                rho_set=rho_set, &
     884              :                                deriv_set=deriv_set, &
     885       133939 :                                deriv_order=deriv_order)
     886              : 
     887       133939 :       CALL divide_by_norm_drho(deriv_set, rho_set, lsd)
     888              : 
     889       133939 :       CALL timestop(handle)
     890              : 
     891       133939 :    END SUBROUTINE xc_rho_set_and_dset_create
     892              : 
     893              : ! **************************************************************************************************
     894              : !> \brief smooths the cutoff on rho with a function smoothderiv_rho that is 0
     895              : !>      for rho<rho_cutoff and 1 for rho>rho_cutoff*rho_smooth_cutoff_range:
     896              : !>      E= integral e_0*smoothderiv_rho => dE/d...= de/d... * smooth,
     897              : !>      dE/drho = de/drho * smooth + e_0 * dsmooth/drho
     898              : !> \param pot the potential to smooth
     899              : !> \param rho , rhoa,rhob: the value of the density (used to apply the cutoff)
     900              : !> \param rhoa ...
     901              : !> \param rhob ...
     902              : !> \param rho_cutoff the value at whch the cutoff function must go to 0
     903              : !> \param rho_smooth_cutoff_range range of the smoothing
     904              : !> \param e_0 value of e_0, if given it is assumed that pot is the derivative
     905              : !>        wrt. to rho, and needs the dsmooth*e_0 contribution
     906              : !> \param e_0_scale_factor ...
     907              : !> \author Fawzi Mohamed
     908              : ! **************************************************************************************************
     909       254857 :    SUBROUTINE smooth_cutoff(pot, rho, rhoa, rhob, rho_cutoff, &
     910              :                             rho_smooth_cutoff_range, e_0, e_0_scale_factor)
     911              :       REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
     912              :          POINTER                                         :: pot, rho, rhoa, rhob
     913              :       REAL(kind=dp), INTENT(in)                          :: rho_cutoff, rho_smooth_cutoff_range
     914              :       REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, &
     915              :          POINTER                                         :: e_0
     916              :       REAL(kind=dp), INTENT(in), OPTIONAL                :: e_0_scale_factor
     917              : 
     918              :       INTEGER                                            :: i, j, k
     919              :       INTEGER, DIMENSION(2, 3)                           :: bo
     920              :       REAL(kind=dp) :: my_e_0_scale_factor, my_rho, my_rho_n, my_rho_n2, rho_smooth_cutoff, &
     921              :                        rho_smooth_cutoff_2, rho_smooth_cutoff_range_2
     922              : 
     923       254857 :       CPASSERT(ASSOCIATED(pot))
     924      1019428 :       bo(1, :) = LBOUND(pot)
     925      1019428 :       bo(2, :) = UBOUND(pot)
     926       254857 :       my_e_0_scale_factor = 1.0_dp
     927       254857 :       IF (PRESENT(e_0_scale_factor)) my_e_0_scale_factor = e_0_scale_factor
     928       254857 :       rho_smooth_cutoff = rho_cutoff*rho_smooth_cutoff_range
     929       254857 :       rho_smooth_cutoff_2 = (rho_cutoff + rho_smooth_cutoff)/2
     930       254857 :       rho_smooth_cutoff_range_2 = rho_smooth_cutoff_2 - rho_cutoff
     931              : 
     932       254857 :       IF (rho_smooth_cutoff_range > 0.0_dp) THEN
     933            2 :          IF (PRESENT(e_0)) THEN
     934            0 :             CPASSERT(ASSOCIATED(e_0))
     935            0 :             IF (ASSOCIATED(rho)) THEN
     936              : !$OMP PARALLEL DO DEFAULT(NONE) &
     937              : !$OMP             SHARED(bo,e_0,pot,rho,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
     938              : !$OMP                    rho_smooth_cutoff_range_2,my_e_0_scale_factor) &
     939              : !$OMP             PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
     940            0 : !$OMP             COLLAPSE(3)
     941              :                DO k = bo(1, 3), bo(2, 3)
     942              :                   DO j = bo(1, 2), bo(2, 2)
     943              :                      DO i = bo(1, 1), bo(2, 1)
     944              :                         my_rho = rho(i, j, k)
     945              :                         IF (my_rho < rho_smooth_cutoff) THEN
     946              :                            IF (my_rho < rho_cutoff) THEN
     947              :                               pot(i, j, k) = 0.0_dp
     948              :                            ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
     949              :                               my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
     950              :                               my_rho_n2 = my_rho_n*my_rho_n
     951              :                               pot(i, j, k) = pot(i, j, k)* &
     952              :                                              my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
     953              :                                              my_e_0_scale_factor*e_0(i, j, k)* &
     954              :                                              my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
     955              :                                              /rho_smooth_cutoff_range_2
     956              :                            ELSE
     957              :                               my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
     958              :                               my_rho_n2 = my_rho_n*my_rho_n
     959              :                               pot(i, j, k) = pot(i, j, k)* &
     960              :                                              (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
     961              :                                              + my_e_0_scale_factor*e_0(i, j, k)* &
     962              :                                              my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
     963              :                                              /rho_smooth_cutoff_range_2
     964              :                            END IF
     965              :                         END IF
     966              :                      END DO
     967              :                   END DO
     968              :                END DO
     969              : !$OMP END PARALLEL DO
     970              :             ELSE
     971              : !$OMP PARALLEL DO DEFAULT(NONE) &
     972              : !$OMP             SHARED(bo,pot,e_0,rhoa,rhob,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
     973              : !$OMP                    rho_smooth_cutoff_range_2,my_e_0_scale_factor) &
     974              : !$OMP             PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
     975            0 : !$OMP             COLLAPSE(3)
     976              :                DO k = bo(1, 3), bo(2, 3)
     977              :                   DO j = bo(1, 2), bo(2, 2)
     978              :                      DO i = bo(1, 1), bo(2, 1)
     979              :                         my_rho = rhoa(i, j, k) + rhob(i, j, k)
     980              :                         IF (my_rho < rho_smooth_cutoff) THEN
     981              :                            IF (my_rho < rho_cutoff) THEN
     982              :                               pot(i, j, k) = 0.0_dp
     983              :                            ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
     984              :                               my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
     985              :                               my_rho_n2 = my_rho_n*my_rho_n
     986              :                               pot(i, j, k) = pot(i, j, k)* &
     987              :                                              my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2) + &
     988              :                                              my_e_0_scale_factor*e_0(i, j, k)* &
     989              :                                              my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
     990              :                                              /rho_smooth_cutoff_range_2
     991              :                            ELSE
     992              :                               my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
     993              :                               my_rho_n2 = my_rho_n*my_rho_n
     994              :                               pot(i, j, k) = pot(i, j, k)* &
     995              :                                              (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)) &
     996              :                                              + my_e_0_scale_factor*e_0(i, j, k)* &
     997              :                                              my_rho_n2*(3.0_dp - 2.0_dp*my_rho_n) &
     998              :                                              /rho_smooth_cutoff_range_2
     999              :                            END IF
    1000              :                         END IF
    1001              :                      END DO
    1002              :                   END DO
    1003              :                END DO
    1004              : !$OMP END PARALLEL DO
    1005              :             END IF
    1006              :          ELSE
    1007            2 :             IF (ASSOCIATED(rho)) THEN
    1008              : !$OMP PARALLEL DO DEFAULT(NONE) &
    1009              : !$OMP             SHARED(bo,pot,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
    1010              : !$OMP                    rho_smooth_cutoff_range_2,rho) &
    1011              : !$OMP             PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
    1012            2 : !$OMP             COLLAPSE(3)
    1013              :                DO k = bo(1, 3), bo(2, 3)
    1014              :                   DO j = bo(1, 2), bo(2, 2)
    1015              :                      DO i = bo(1, 1), bo(2, 1)
    1016              :                         my_rho = rho(i, j, k)
    1017              :                         IF (my_rho < rho_smooth_cutoff) THEN
    1018              :                            IF (my_rho < rho_cutoff) THEN
    1019              :                               pot(i, j, k) = 0.0_dp
    1020              :                            ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
    1021              :                               my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
    1022              :                               my_rho_n2 = my_rho_n*my_rho_n
    1023              :                               pot(i, j, k) = pot(i, j, k)* &
    1024              :                                              my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
    1025              :                            ELSE
    1026              :                               my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
    1027              :                               my_rho_n2 = my_rho_n*my_rho_n
    1028              :                               pot(i, j, k) = pot(i, j, k)* &
    1029              :                                              (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
    1030              :                            END IF
    1031              :                         END IF
    1032              :                      END DO
    1033              :                   END DO
    1034              :                END DO
    1035              : !$OMP END PARALLEL DO
    1036              :             ELSE
    1037              : !$OMP PARALLEL DO DEFAULT(NONE) &
    1038              : !$OMP             SHARED(bo,pot,rho_cutoff,rho_smooth_cutoff,rho_smooth_cutoff_2, &
    1039              : !$OMP                    rho_smooth_cutoff_range_2,rhoa,rhob) &
    1040              : !$OMP             PRIVATE(k,j,i,my_rho,my_rho_n,my_rho_n2) &
    1041            0 : !$OMP             COLLAPSE(3)
    1042              :                DO k = bo(1, 3), bo(2, 3)
    1043              :                   DO j = bo(1, 2), bo(2, 2)
    1044              :                      DO i = bo(1, 1), bo(2, 1)
    1045              :                         my_rho = rhoa(i, j, k) + rhob(i, j, k)
    1046              :                         IF (my_rho < rho_smooth_cutoff) THEN
    1047              :                            IF (my_rho < rho_cutoff) THEN
    1048              :                               pot(i, j, k) = 0.0_dp
    1049              :                            ELSEIF (my_rho < rho_smooth_cutoff_2) THEN
    1050              :                               my_rho_n = (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
    1051              :                               my_rho_n2 = my_rho_n*my_rho_n
    1052              :                               pot(i, j, k) = pot(i, j, k)* &
    1053              :                                              my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2)
    1054              :                            ELSE
    1055              :                               my_rho_n = 2.0_dp - (my_rho - rho_cutoff)/rho_smooth_cutoff_range_2
    1056              :                               my_rho_n2 = my_rho_n*my_rho_n
    1057              :                               pot(i, j, k) = pot(i, j, k)* &
    1058              :                                              (1.0_dp - my_rho_n2*(my_rho_n - 0.5_dp*my_rho_n2))
    1059              :                            END IF
    1060              :                         END IF
    1061              :                      END DO
    1062              :                   END DO
    1063              :                END DO
    1064              : !$OMP END PARALLEL DO
    1065              :             END IF
    1066              :          END IF
    1067              :       END IF
    1068       254857 :    END SUBROUTINE smooth_cutoff
    1069              : 
    1070           64 :    SUBROUTINE calc_xc_density(pot, rho, rho_cutoff)
    1071              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: pot
    1072              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT)         :: rho
    1073              :       REAL(kind=dp), INTENT(in)                          :: rho_cutoff
    1074              : 
    1075              :       INTEGER                                            :: i, j, k, nspins
    1076              :       INTEGER, DIMENSION(2, 3)                           :: bo
    1077              :       REAL(kind=dp)                                      :: eps1, eps2, my_rho, my_pot
    1078              : 
    1079          256 :       bo(1, :) = LBOUND(pot%array)
    1080          256 :       bo(2, :) = UBOUND(pot%array)
    1081           64 :       nspins = SIZE(rho)
    1082              : 
    1083           64 :       eps1 = rho_cutoff*1.E-4_dp
    1084           64 :       eps2 = rho_cutoff
    1085              : 
    1086         3160 :       DO k = bo(1, 3), bo(2, 3)
    1087       161272 :          DO j = bo(1, 2), bo(2, 2)
    1088      4529376 :             DO i = bo(1, 1), bo(2, 1)
    1089      4368168 :                my_pot = pot%array(i, j, k)
    1090      4368168 :                IF (nspins == 2) THEN
    1091       339714 :                   my_rho = rho(1)%array(i, j, k) + rho(2)%array(i, j, k)
    1092              :                ELSE
    1093      4028454 :                   my_rho = rho(1)%array(i, j, k)
    1094              :                END IF
    1095      4526280 :                IF (my_rho > eps1) THEN
    1096      4139578 :                   pot%array(i, j, k) = my_pot/my_rho
    1097       228590 :                ELSE IF (my_rho < eps2) THEN
    1098       228590 :                   pot%array(i, j, k) = 0.0_dp
    1099              :                ELSE
    1100            0 :                   pot%array(i, j, k) = MIN(my_pot/my_rho, my_rho**(1._dp/3._dp))
    1101              :                END IF
    1102              :             END DO
    1103              :          END DO
    1104              :       END DO
    1105              : 
    1106           64 :    END SUBROUTINE calc_xc_density
    1107              : 
    1108              : ! **************************************************************************************************
    1109              : !> \brief Exchange and Correlation functional calculations
    1110              : !> \param vxc_rho will contain the v_xc part that depend on rho
    1111              : !>        (if one of the chosen xc functionals has it it is allocated and you
    1112              : !>        are responsible for it)
    1113              : !> \param vxc_tau will contain the kinetic tau part of v_xc
    1114              : !>        (if one of the chosen xc functionals has it it is allocated and you
    1115              : !>        are responsible for it)
    1116              : !> \param exc the xc energy
    1117              : !> \param rho_r the value of the density in the real space
    1118              : !> \param rho_g value of the density in the g space (needs to be associated
    1119              : !>        only for gradient corrections)
    1120              : !> \param tau value of the kinetic density tau on the grid (can be null,
    1121              : !>        used only with meta functionals)
    1122              : !> \param xc_section which functional to calculate, and how to do it
    1123              : !> \param pw_pool the pool for the grids
    1124              : !> \param compute_virial ...
    1125              : !> \param virial_xc ...
    1126              : !> \param exc_r the value of the xc functional in the real space
    1127              : !> \par History
    1128              : !>      JGH (13-Jun-2002): adaptation to new functionals
    1129              : !>      Fawzi (11.2002): drho_g(1:3)->drho_g
    1130              : !>      Fawzi (1.2003). lsd version
    1131              : !>      Fawzi (11.2003): version using the new xc interface
    1132              : !>      Fawzi (03.2004): fft free for smoothed density and derivs, gga lsd
    1133              : !>      Fawzi (04.2004): metafunctionals
    1134              : !>      mguidon (12.2008) : laplace functionals
    1135              : !> \author fawzi; based LDA version of JGH, based on earlier version of apsi
    1136              : !> \note
    1137              : !>      Beware: some really dirty pointer handling!
    1138              : !>      energy should be kept consistent with xc_exc_calc
    1139              : ! **************************************************************************************************
    1140       111239 :    SUBROUTINE xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, &
    1141              :                                pw_pool, compute_virial, virial_xc, exc_r)
    1142              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: vxc_rho, vxc_tau
    1143              :       REAL(KIND=dp), INTENT(out)                         :: exc
    1144              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               :: rho_r, tau
    1145              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
    1146              :       TYPE(section_vals_type), POINTER                   :: xc_section
    1147              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
    1148              :       LOGICAL                                            :: compute_virial
    1149              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: virial_xc
    1150              :       TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL            :: exc_r
    1151              : 
    1152              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_vxc_pw_create'
    1153              :       INTEGER, DIMENSION(2), PARAMETER :: norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob]
    1154              : 
    1155              :       INTEGER                                            :: handle, idir, ispin, jdir, &
    1156              :                                                             npoints, nspins, &
    1157              :                                                             xc_deriv_method_id, xc_rho_smooth_id, deriv_id
    1158              :       INTEGER, DIMENSION(2, 3)                           :: bo
    1159              :       LOGICAL                                            :: dealloc_pw_to_deriv, has_laplace, &
    1160              :                                                             has_tau, lsd, use_virial, has_gradient, &
    1161              :                                                             has_derivs, has_rho, dealloc_pw_to_deriv_rho
    1162              :       REAL(KIND=dp)                                      :: density_smooth_cut_range, drho_cutoff, &
    1163              :                                                             rho_cutoff
    1164       111239 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: deriv_data, norm_drho, norm_drho_spin, &
    1165       222478 :                                                             rho, rhoa, rhob
    1166              :       TYPE(cp_sll_xc_deriv_type), POINTER                :: pos
    1167              :       TYPE(pw_grid_type), POINTER                        :: pw_grid
    1168       778673 :       TYPE(pw_r3d_rs_type), DIMENSION(3)                      :: pw_to_deriv, pw_to_deriv_rho
    1169              :       TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
    1170              :       TYPE(pw_r3d_rs_type)                                      ::  v_drho_r, virial_pw
    1171              :       TYPE(xc_derivative_set_type)                       :: deriv_set
    1172              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att
    1173              :       TYPE(xc_rho_set_type)                              :: rho_set
    1174              : 
    1175       111239 :       CALL timeset(routineN, handle)
    1176       111239 :       NULLIFY (norm_drho_spin, norm_drho, pos)
    1177              : 
    1178       111239 :       pw_grid => rho_r(1)%pw_grid
    1179              : 
    1180       111239 :       CPASSERT(ASSOCIATED(xc_section))
    1181       111239 :       CPASSERT(ASSOCIATED(pw_pool))
    1182       111239 :       CPASSERT(.NOT. ASSOCIATED(vxc_rho))
    1183       111239 :       CPASSERT(.NOT. ASSOCIATED(vxc_tau))
    1184       111239 :       nspins = SIZE(rho_r)
    1185       111239 :       lsd = (nspins /= 1)
    1186       111239 :       IF (lsd) THEN
    1187        22129 :          CPASSERT(nspins == 2)
    1188              :       END IF
    1189              : 
    1190       111239 :       use_virial = compute_virial
    1191       111239 :       virial_xc = 0.0_dp
    1192              : 
    1193      1112390 :       bo = rho_r(1)%pw_grid%bounds_local
    1194       111239 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
    1195              : 
    1196              :       ! calculate the potential derivatives
    1197              :       CALL xc_rho_set_and_dset_create(rho_set=rho_set, deriv_set=deriv_set, &
    1198              :                                       deriv_order=1, rho_r=rho_r, rho_g=rho_g, tau=tau, &
    1199              :                                       xc_section=xc_section, &
    1200              :                                       pw_pool=pw_pool, &
    1201       111239 :                                       calc_potential=.TRUE.)
    1202              : 
    1203              :       CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
    1204       111239 :                                 i_val=xc_deriv_method_id)
    1205              :       CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
    1206       111239 :                                 i_val=xc_rho_smooth_id)
    1207              :       CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
    1208       111239 :                                 r_val=density_smooth_cut_range)
    1209              : 
    1210              :       CALL xc_rho_set_get(rho_set, rho_cutoff=rho_cutoff, &
    1211       111239 :                           drho_cutoff=drho_cutoff)
    1212              : 
    1213       111239 :       CALL check_for_derivatives(deriv_set, lsd, has_rho, has_gradient, has_tau, has_laplace)
    1214              :       ! check for unknown derivatives
    1215       111239 :       has_derivs = has_rho .OR. has_gradient .OR. has_tau .OR. has_laplace
    1216              : 
    1217       467085 :       ALLOCATE (vxc_rho(nspins))
    1218              : 
    1219              :       CALL xc_rho_set_get(rho_set, rho=rho, rhoa=rhoa, rhob=rhob, &
    1220       111239 :                           can_return_null=.TRUE.)
    1221              : 
    1222              :       ! recover the vxc arrays
    1223       111239 :       IF (lsd) THEN
    1224        22129 :          CALL xc_dset_recover_pw(deriv_set, [deriv_rhoa], vxc_rho(1), pw_grid, pw_pool)
    1225        22129 :          CALL xc_dset_recover_pw(deriv_set, [deriv_rhob], vxc_rho(2), pw_grid, pw_pool)
    1226              : 
    1227              :       ELSE
    1228        89110 :          CALL xc_dset_recover_pw(deriv_set, [deriv_rho], vxc_rho(1), pw_grid, pw_pool)
    1229              :       END IF
    1230              : 
    1231       111239 :       deriv_att => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
    1232       111239 :       IF (ASSOCIATED(deriv_att)) THEN
    1233        60013 :          CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
    1234              : 
    1235              :          CALL xc_rho_set_get(rho_set, norm_drho=norm_drho, &
    1236              :                              rho_cutoff=rho_cutoff, &
    1237              :                              drho_cutoff=drho_cutoff, &
    1238        60013 :                              can_return_null=.TRUE.)
    1239        60013 :          CALL xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv_rho, drho=pw_to_deriv_rho)
    1240              : 
    1241        60013 :          CPASSERT(ASSOCIATED(deriv_data))
    1242        60013 :          IF (use_virial) THEN
    1243         1568 :             CALL pw_pool%create_pw(virial_pw)
    1244         1568 :             CALL pw_zero(virial_pw)
    1245         6272 :             DO idir = 1, 3
    1246         4704 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(virial_pw,pw_to_deriv_rho,deriv_data,idir)
    1247              :                virial_pw%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
    1248              : !$OMP END PARALLEL WORKSHARE
    1249        15680 :                DO jdir = 1, idir
    1250              :                   virial_xc(idir, jdir) = -pw_grid%dvol* &
    1251              :                                           accurate_dot_product(virial_pw%array(:, :, :), &
    1252         9408 :                                                                pw_to_deriv_rho(jdir)%array(:, :, :))
    1253        14112 :                   virial_xc(jdir, idir) = virial_xc(idir, jdir)
    1254              :                END DO
    1255              :             END DO
    1256         1568 :             CALL pw_pool%give_back_pw(virial_pw)
    1257              :          END IF ! use_virial
    1258       240052 :          DO idir = 1, 3
    1259       180039 :             CPASSERT(ASSOCIATED(pw_to_deriv_rho(idir)%array))
    1260       240052 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv_rho,idir)
    1261              :             pw_to_deriv_rho(idir)%array(:, :, :) = pw_to_deriv_rho(idir)%array(:, :, :)*deriv_data(:, :, :)
    1262              : !$OMP END PARALLEL WORKSHARE
    1263              :          END DO
    1264              : 
    1265              :          ! Deallocate pw to save memory
    1266        60013 :          CALL pw_pool%give_back_cr3d(deriv_att%deriv_data)
    1267              : 
    1268              :       END IF
    1269              : 
    1270       111239 :       IF ((has_gradient .AND. xc_requires_tmp_g(xc_deriv_method_id)) .OR. pw_grid%spherical) THEN
    1271        59783 :          CALL pw_pool%create_pw(vxc_g)
    1272        59783 :          IF (.NOT. pw_grid%spherical) THEN
    1273        59783 :             CALL pw_pool%create_pw(tmp_g)
    1274              :          END IF
    1275              :       END IF
    1276              : 
    1277       244607 :       DO ispin = 1, nspins
    1278              : 
    1279       133368 :          IF (lsd) THEN
    1280        44258 :             IF (ispin == 1) THEN
    1281              :                CALL xc_rho_set_get(rho_set, norm_drhoa=norm_drho_spin, &
    1282        22129 :                                    can_return_null=.TRUE.)
    1283        22129 :                IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
    1284        14080 :                   rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhoa=pw_to_deriv)
    1285              :             ELSE
    1286              :                CALL xc_rho_set_get(rho_set, norm_drhob=norm_drho_spin, &
    1287        22129 :                                    can_return_null=.TRUE.)
    1288        22129 :                IF (ASSOCIATED(norm_drho_spin)) CALL xc_rho_set_recover_pw( &
    1289        14080 :                   rho_set, pw_grid, pw_pool, dealloc_pw_to_deriv, drhob=pw_to_deriv)
    1290              :             END IF
    1291              : 
    1292        88516 :             deriv_att => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)])
    1293        44258 :             IF (ASSOCIATED(deriv_att)) THEN
    1294              :                CPASSERT(lsd)
    1295        28160 :                CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
    1296              : 
    1297        28160 :                IF (use_virial) THEN
    1298          104 :                   CALL pw_pool%create_pw(virial_pw)
    1299          104 :                   CALL pw_zero(virial_pw)
    1300          416 :                   DO idir = 1, 3
    1301          312 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,pw_to_deriv,virial_pw,idir)
    1302              :                      virial_pw%array(:, :, :) = pw_to_deriv(idir)%array(:, :, :)*deriv_data(:, :, :)
    1303              : !$OMP END PARALLEL WORKSHARE
    1304         1040 :                      DO jdir = 1, idir
    1305              :                         virial_xc(idir, jdir) = virial_xc(idir, jdir) - pw_grid%dvol* &
    1306              :                                                 accurate_dot_product(virial_pw%array(:, :, :), &
    1307          624 :                                                                      pw_to_deriv(jdir)%array(:, :, :))
    1308          936 :                         virial_xc(jdir, idir) = virial_xc(idir, jdir)
    1309              :                      END DO
    1310              :                   END DO
    1311          104 :                   CALL pw_pool%give_back_pw(virial_pw)
    1312              :                END IF ! use_virial
    1313              : 
    1314       112640 :                DO idir = 1, 3
    1315       112640 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_data,idir,pw_to_deriv)
    1316              :                   pw_to_deriv(idir)%array(:, :, :) = deriv_data(:, :, :)*pw_to_deriv(idir)%array(:, :, :)
    1317              : !$OMP END PARALLEL WORKSHARE
    1318              :                END DO
    1319              :             END IF ! deriv_att
    1320              : 
    1321              :          END IF ! LSD
    1322              : 
    1323       133368 :          IF (ASSOCIATED(pw_to_deriv_rho(1)%array)) THEN
    1324        73335 :             IF (.NOT. ASSOCIATED(pw_to_deriv(1)%array)) THEN
    1325        46691 :                pw_to_deriv = pw_to_deriv_rho
    1326        46691 :                dealloc_pw_to_deriv = ((.NOT. lsd) .OR. (ispin == 2))
    1327        46691 :                dealloc_pw_to_deriv = dealloc_pw_to_deriv .AND. dealloc_pw_to_deriv_rho
    1328              :             ELSE
    1329              :                ! This branch is called in case of open-shell systems
    1330              :                ! Add the contributions from norm_drho and norm_drho_spin
    1331       106576 :                DO idir = 1, 3
    1332        79932 :                   CALL pw_axpy(pw_to_deriv_rho(idir), pw_to_deriv(idir))
    1333       106576 :                   IF (ispin == 2) THEN
    1334        39966 :                      IF (dealloc_pw_to_deriv_rho) THEN
    1335        39966 :                         CALL pw_pool%give_back_pw(pw_to_deriv_rho(idir))
    1336              :                      END IF
    1337              :                   END IF
    1338              :                END DO
    1339              :             END IF
    1340              :          END IF
    1341              : 
    1342       133368 :          IF (ASSOCIATED(pw_to_deriv(1)%array)) THEN
    1343       299404 :             DO idir = 1, 3
    1344       299404 :                CALL pw_scale(pw_to_deriv(idir), -1.0_dp)
    1345              :             END DO
    1346              : 
    1347        74851 :             CALL xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_rho(ispin))
    1348              : 
    1349        74851 :             IF (dealloc_pw_to_deriv) THEN
    1350       299404 :             DO idir = 1, 3
    1351       299404 :                CALL pw_pool%give_back_pw(pw_to_deriv(idir))
    1352              :             END DO
    1353              :             END IF
    1354              :          END IF
    1355              : 
    1356              :          ! Add laplace part to vxc_rho
    1357       133368 :          IF (has_laplace) THEN
    1358         1092 :             IF (lsd) THEN
    1359          472 :                IF (ispin == 1) THEN
    1360              :                   deriv_id = deriv_laplace_rhoa
    1361              :                ELSE
    1362          236 :                   deriv_id = deriv_laplace_rhob
    1363              :                END IF
    1364              :             ELSE
    1365              :                deriv_id = deriv_laplace_rho
    1366              :             END IF
    1367              : 
    1368         2184 :             CALL xc_dset_recover_pw(deriv_set, [deriv_id], pw_to_deriv(1), pw_grid)
    1369              : 
    1370         1092 :             IF (use_virial) CALL virial_laplace(rho_r(ispin), pw_pool, virial_xc, pw_to_deriv(1)%array)
    1371              : 
    1372         1092 :             CALL xc_pw_laplace(pw_to_deriv(1), pw_pool, xc_deriv_method_id)
    1373              : 
    1374         1092 :             CALL pw_axpy(pw_to_deriv(1), vxc_rho(ispin))
    1375              : 
    1376         1092 :             CALL pw_pool%give_back_pw(pw_to_deriv(1))
    1377              :          END IF
    1378              : 
    1379       133368 :          IF (pw_grid%spherical) THEN
    1380              :             ! filter vxc
    1381            0 :             CALL pw_transfer(vxc_rho(ispin), vxc_g)
    1382            0 :             CALL pw_transfer(vxc_g, vxc_rho(ispin))
    1383              :          END IF
    1384              :          CALL smooth_cutoff(pot=vxc_rho(ispin)%array, rho=rho, rhoa=rhoa, rhob=rhob, &
    1385              :                             rho_cutoff=rho_cutoff*density_smooth_cut_range, &
    1386       133368 :                             rho_smooth_cutoff_range=density_smooth_cut_range)
    1387              : 
    1388       133368 :          v_drho_r = vxc_rho(ispin)
    1389       133368 :          CALL pw_pool%create_pw(vxc_rho(ispin))
    1390       133368 :          CALL xc_pw_smooth(v_drho_r, vxc_rho(ispin), xc_rho_smooth_id)
    1391       244607 :          CALL pw_pool%give_back_pw(v_drho_r)
    1392              :       END DO
    1393              : 
    1394       111239 :       CALL pw_pool%give_back_pw(vxc_g)
    1395       111239 :       CALL pw_pool%give_back_pw(tmp_g)
    1396              : 
    1397              :       ! 0-deriv -> value of exc
    1398              :       ! this has to be kept consistent with xc_exc_calc
    1399       111239 :       IF (has_derivs) THEN
    1400       110921 :          CALL xc_dset_recover_pw(deriv_set, [INTEGER::], v_drho_r, pw_grid)
    1401              : 
    1402              :          CALL smooth_cutoff(pot=v_drho_r%array, rho=rho, rhoa=rhoa, rhob=rhob, &
    1403              :                             rho_cutoff=rho_cutoff, &
    1404       110921 :                             rho_smooth_cutoff_range=density_smooth_cut_range)
    1405              : 
    1406       110921 :          exc = pw_integrate_function(v_drho_r)
    1407              :          !
    1408              :          ! return the xc functional value at the grid points
    1409              :          !
    1410       110921 :          IF (PRESENT(exc_r)) THEN
    1411           96 :             exc_r = v_drho_r
    1412              :          ELSE
    1413       110825 :             CALL v_drho_r%release()
    1414              :          END IF
    1415              :       ELSE
    1416          318 :          exc = 0.0_dp
    1417              :       END IF
    1418              : 
    1419       111239 :       CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
    1420              : 
    1421              :       ! tau part
    1422       111239 :       IF (has_tau) THEN
    1423         8230 :          ALLOCATE (vxc_tau(nspins))
    1424         2532 :          IF (lsd) THEN
    1425          634 :             CALL xc_dset_recover_pw(deriv_set, [deriv_tau_a], vxc_tau(1), pw_grid)
    1426          634 :             CALL xc_dset_recover_pw(deriv_set, [deriv_tau_b], vxc_tau(2), pw_grid)
    1427              :          ELSE
    1428         1898 :             CALL xc_dset_recover_pw(deriv_set, [deriv_tau], vxc_tau(1), pw_grid)
    1429              :          END IF
    1430         5698 :          DO ispin = 1, nspins
    1431         5698 :             CPASSERT(ASSOCIATED(vxc_tau(ispin)%array))
    1432              :          END DO
    1433              :       END IF
    1434       111239 :       CALL xc_dset_release(deriv_set)
    1435              : 
    1436       111239 :       CALL timestop(handle)
    1437              : 
    1438      2336019 :    END SUBROUTINE xc_vxc_pw_create
    1439              : 
    1440              : ! **************************************************************************************************
    1441              : !> \brief calculates just the exchange and correlation energy
    1442              : !>      (no vxc)
    1443              : !> \param rho_r      realspace density on the grid
    1444              : !> \param rho_g      g-space density on the grid
    1445              : !> \param tau        kinetic energy density on the grid
    1446              : !> \param xc_section XC parameters
    1447              : !> \param pw_pool    pool of plain-wave grids
    1448              : !> \return the XC energy
    1449              : !> \par History
    1450              : !>      11.2003 created [fawzi]
    1451              : !> \author fawzi
    1452              : !> \note
    1453              : !>      has to be kept consistent with xc_vxc_pw_create
    1454              : ! **************************************************************************************************
    1455        21132 :    FUNCTION xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool) &
    1456              :       RESULT(exc)
    1457              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               :: rho_r, tau
    1458              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
    1459              :       TYPE(section_vals_type), POINTER                   :: xc_section
    1460              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
    1461              :       REAL(kind=dp)                                      :: exc
    1462              : 
    1463              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_exc_calc'
    1464              : 
    1465              :       INTEGER                                            :: handle
    1466              :       REAL(dp)                                           :: density_smooth_cut_range, rho_cutoff
    1467        10566 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: e_0
    1468              :       TYPE(xc_derivative_set_type)                       :: deriv_set
    1469              :       TYPE(xc_derivative_type), POINTER                  :: deriv
    1470              :       TYPE(xc_rho_set_type)                              :: rho_set
    1471              : 
    1472        10566 :       CALL timeset(routineN, handle)
    1473        10566 :       NULLIFY (deriv, e_0)
    1474        10566 :       exc = 0.0_dp
    1475              : 
    1476              :       ! this has to be consistent with what is done in xc_vxc_pw_create
    1477              :       CALL xc_rho_set_and_dset_create(rho_set=rho_set, &
    1478              :                                       deriv_set=deriv_set, deriv_order=0, &
    1479              :                                       rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
    1480              :                                       pw_pool=pw_pool, &
    1481        10566 :                                       calc_potential=.FALSE.)
    1482        10566 :       deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
    1483              : 
    1484        10566 :       IF (ASSOCIATED(deriv)) THEN
    1485        10566 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
    1486              : 
    1487              :          CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
    1488        10566 :                                    r_val=rho_cutoff)
    1489              :          CALL section_vals_val_get(xc_section, "DENSITY_SMOOTH_CUTOFF_RANGE", &
    1490        10566 :                                    r_val=density_smooth_cut_range)
    1491              :          CALL smooth_cutoff(pot=e_0, rho=rho_set%rho, &
    1492              :                             rhoa=rho_set%rhoa, rhob=rho_set%rhob, &
    1493              :                             rho_cutoff=rho_cutoff, &
    1494        10566 :                             rho_smooth_cutoff_range=density_smooth_cut_range)
    1495              : 
    1496        10566 :          exc = accurate_sum(e_0)*rho_r(1)%pw_grid%dvol
    1497        10566 :          IF (rho_r(1)%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN
    1498        10428 :             CALL rho_r(1)%pw_grid%para%group%sum(exc)
    1499              :          END IF
    1500              : 
    1501        10566 :          CALL xc_rho_set_release(rho_set, pw_pool=pw_pool)
    1502        10566 :          CALL xc_dset_release(deriv_set)
    1503              :       END IF
    1504        10566 :       CALL timestop(handle)
    1505       200754 :    END FUNCTION xc_exc_calc
    1506              : 
    1507              : ! **************************************************************************************************
    1508              : !> \brief Caller routine to calculate the second order potential in the direction of rho1_r
    1509              : !> \param v_xc XC potential, will be allocated, to be integrated with the KS density
    1510              : !> \param v_xc_tau ...
    1511              : !> \param deriv_set XC derivatives from xc_prep_2nd_deriv
    1512              : !> \param rho_set XC rho set from KS rho from xc_prep_2nd_deriv
    1513              : !> \param rho1_r first-order density in r space
    1514              : !> \param rho1_g first-order density in g space
    1515              : !> \param tau1_r ...
    1516              : !> \param pw_pool pw pool to create new grids
    1517              : !> \param xc_section XC section to calculate the derivatives from
    1518              : !> \param gapw whether to carry out GAPW (not possible with numerical derivatives)
    1519              : !> \param vxg GAPW potential
    1520              : !> \param do_excitations ...
    1521              : !> \param do_triplet ...
    1522              : !> \param compute_virial ...
    1523              : !> \param virial_xc virial terms will be collected here
    1524              : ! **************************************************************************************************
    1525        13900 :    SUBROUTINE xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, &
    1526              :                                 pw_pool, xc_section, gapw, vxg, &
    1527              :                                 do_excitations, do_triplet, &
    1528              :                                 compute_virial, virial_xc)
    1529              : 
    1530              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_xc, v_xc_tau
    1531              :       TYPE(xc_derivative_set_type)                       :: deriv_set
    1532              :       TYPE(xc_rho_set_type)                              :: rho_set
    1533              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, tau1_r
    1534              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_g
    1535              :       TYPE(pw_pool_type), INTENT(IN), POINTER            :: pw_pool
    1536              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: xc_section
    1537              :       LOGICAL, INTENT(IN)                                :: gapw
    1538              :       REAL(KIND=dp), DIMENSION(:, :, :, :), OPTIONAL, &
    1539              :          POINTER                                         :: vxg
    1540              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_excitations, &
    1541              :                                                             do_triplet, compute_virial
    1542              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
    1543              :          OPTIONAL                                        :: virial_xc
    1544              : 
    1545              :       CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv'
    1546              : 
    1547              :       INTEGER                                            :: handle, ispin, nspins
    1548              :       INTEGER, DIMENSION(2, 3)                           :: bo
    1549              :       LOGICAL                                            :: lsd, my_compute_virial, &
    1550              :                                                             my_do_excitations, &
    1551              :                                                             my_do_triplet
    1552              :       REAL(KIND=dp)                                      :: fac
    1553              :       TYPE(section_vals_type), POINTER                   :: xc_fun_section
    1554              :       TYPE(xc_rho_cflags_type)                           :: needs
    1555              :       TYPE(xc_rho_set_type)                              :: rho1_set
    1556              : 
    1557        13900 :       CALL timeset(routineN, handle)
    1558              : 
    1559        13900 :       my_compute_virial = .FALSE.
    1560        13900 :       IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
    1561              : 
    1562        13900 :       my_do_excitations = .FALSE.
    1563        13900 :       IF (PRESENT(do_excitations)) my_do_excitations = do_excitations
    1564              : 
    1565        13900 :       my_do_triplet = .FALSE.
    1566        13900 :       IF (PRESENT(do_triplet)) my_do_triplet = do_triplet
    1567              : 
    1568        13900 :       nspins = SIZE(rho1_r)
    1569        13900 :       lsd = (nspins == 2)
    1570        13900 :       IF (nspins == 1 .AND. my_do_excitations .AND. my_do_triplet) THEN
    1571            0 :          nspins = 2
    1572            0 :          lsd = .TRUE.
    1573              :       END IF
    1574              : 
    1575        13900 :       NULLIFY (v_xc, v_xc_tau)
    1576        57456 :       ALLOCATE (v_xc(nspins))
    1577        29656 :       DO ispin = 1, nspins
    1578        15756 :          CALL pw_pool%create_pw(v_xc(ispin))
    1579        29656 :          CALL pw_zero(v_xc(ispin))
    1580              :       END DO
    1581              : 
    1582        13900 :       xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
    1583        13900 :       needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
    1584              : 
    1585        13900 :       IF (needs%tau .OR. needs%tau_spin) THEN
    1586          536 :          IF (.NOT. ASSOCIATED(tau1_r)) &
    1587            0 :             CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
    1588         1736 :          ALLOCATE (v_xc_tau(nspins))
    1589         1200 :          DO ispin = 1, nspins
    1590          664 :             CALL pw_pool%create_pw(v_xc_tau(ispin))
    1591        14564 :             CALL pw_zero(v_xc_tau(ispin))
    1592              :          END DO
    1593              :       END IF
    1594              : 
    1595        13900 :       IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
    1596              :          !------!
    1597              :          ! rho1 !
    1598              :          !------!
    1599       136000 :          bo = rho1_r(1)%pw_grid%bounds_local
    1600              :          ! create the place where to store the argument for the functionals
    1601              :          CALL xc_rho_set_create(rho1_set, bo, &
    1602              :                                 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
    1603              :                                 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
    1604        13600 :                                 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
    1605              : 
    1606              :          ! calculate the arguments needed by the functionals
    1607              :          CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
    1608              :                                 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
    1609              :                                 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
    1610        13600 :                                 pw_pool)
    1611              : 
    1612        13600 :          fac = 0._dp
    1613        13600 :          IF (nspins == 1 .AND. my_do_excitations) THEN
    1614         1124 :             IF (my_do_triplet) fac = -1.0_dp
    1615              :          END IF
    1616              : 
    1617              :          CALL xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, &
    1618              :                                            rho1_set, pw_pool, xc_section, gapw, vxg=vxg, &
    1619        13600 :                                            tddfpt_fac=fac, compute_virial=compute_virial, virial_xc=virial_xc)
    1620              : 
    1621        13600 :          CALL xc_rho_set_release(rho1_set)
    1622              : 
    1623              :       ELSE
    1624          300 :          IF (gapw) CPABORT("Numerical 2nd derivatives not implemented with GAPW")
    1625              : 
    1626              :          CALL xc_calc_2nd_deriv_numerical(v_xc, v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
    1627              :                                           pw_pool, xc_section, &
    1628              :                                           my_do_excitations .AND. my_do_triplet, &
    1629          300 :                                           compute_virial, virial_xc, deriv_set)
    1630              :       END IF
    1631              : 
    1632        13900 :       CALL timestop(handle)
    1633              : 
    1634       305800 :    END SUBROUTINE xc_calc_2nd_deriv
    1635              : 
    1636              : ! **************************************************************************************************
    1637              : !> \brief calculates 2nd derivative numerically
    1638              : !> \param v_xc potential to be calculated (has to be allocated already)
    1639              : !> \param v_tau tau-part of the potential to be calculated (has to be allocated already)
    1640              : !> \param rho_set KS density from xc_prep_2nd_deriv
    1641              : !> \param rho1_r first-order density in r-space
    1642              : !> \param rho1_g first-order density in g-space
    1643              : !> \param tau1_r first-order kinetic-energy density in r-space
    1644              : !> \param pw_pool pw pool for new grids
    1645              : !> \param xc_section XC section to calculate the derivatives from
    1646              : !> \param do_triplet ...
    1647              : !> \param calc_virial whether to calculate virial terms
    1648              : !> \param virial_xc collects stress tensor components (no metaGGAs!)
    1649              : !> \param deriv_set deriv set from xc_prep_2nd_deriv (only for virials)
    1650              : ! **************************************************************************************************
    1651          340 :    SUBROUTINE xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, &
    1652              :                                           pw_pool, xc_section, &
    1653              :                                           do_triplet, calc_virial, virial_xc, deriv_set)
    1654              : 
    1655              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER :: v_xc, v_tau
    1656              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
    1657              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER   :: rho1_r, tau1_r
    1658              :       TYPE(pw_c1d_gs_type), DIMENSION(:), INTENT(IN), POINTER :: rho1_g
    1659              :       TYPE(pw_pool_type), INTENT(IN), POINTER            :: pw_pool
    1660              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: xc_section
    1661              :       LOGICAL, INTENT(IN)                                :: do_triplet
    1662              :       LOGICAL, INTENT(IN), OPTIONAL                      :: calc_virial
    1663              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
    1664              :          OPTIONAL                                        :: virial_xc
    1665              :       TYPE(xc_derivative_set_type), OPTIONAL             :: deriv_set
    1666              : 
    1667              :       CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_numerical'
    1668              :       REAL(KIND=dp), DIMENSION(-4:4, 4), PARAMETER :: &
    1669              :          weights = RESHAPE([0.0_dp, 0.0_dp, 0.0_dp, -0.5_dp, 0.0_dp, 0.5_dp, 0.0_dp, 0.0_dp, 0.0_dp, &
    1670              :                            0.0_dp, 0.0_dp, 1.0_dp/12.0_dp, -2.0_dp/3.0_dp, 0.0_dp, 2.0_dp/3.0_dp, -1.0_dp/12.0_dp, 0.0_dp, 0.0_dp, &
    1671              :                             0.0_dp, -1.0_dp/60.0_dp, 0.15_dp, -0.75_dp, 0.0_dp, 0.75_dp, -0.15_dp, 1.0_dp/60.0_dp, 0.0_dp, &
    1672              :             1.0_dp/280.0_dp, -4.0_dp/105.0_dp, 0.2_dp, -0.8_dp, 0.0_dp, 0.8_dp, -0.2_dp, 4.0_dp/105.0_dp, -1.0_dp/280.0_dp], [9, 4])
    1673              : 
    1674              :       INTEGER                                            :: handle, idir, ispin, nspins, istep, nsteps
    1675              :       INTEGER, DIMENSION(2, 3)                           :: bo
    1676              :       LOGICAL                                            :: gradient_f, lsd, my_calc_virial, tau_f, laplace_f, rho_f
    1677              :       REAL(KIND=dp)                                      :: exc, gradient_cut, h, weight, step, rho_cutoff
    1678          340 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dr1dr, dra1dra, drb1drb
    1679              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_dummy
    1680          340 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: norm_drho, norm_drho2, norm_drho2a, &
    1681          340 :                                                             norm_drho2b, norm_drhoa, norm_drhob, &
    1682         1020 :                                                             rho, rho1, rho1a, rho1b, rhoa, rhob, &
    1683          680 :                                                             tau_a, tau_b, tau, tau1, tau1a, tau1b, laplace, laplace1, &
    1684          340 :                                                             laplacea, laplaceb, laplace1a, laplace1b, &
    1685          680 :                                                             laplace2, laplace2a, laplace2b, deriv_data
    1686         8160 :       TYPE(cp_3d_r_cp_type), DIMENSION(3)                :: drho, drho1, drho1a, drho1b, drhoa, drhob
    1687              :       TYPE(pw_r3d_rs_type)                                      :: v_drho, v_drhoa, v_drhob
    1688          340 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: vxc_rho, vxc_tau
    1689          340 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
    1690          340 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               ::  rho_r, tau_r
    1691              :       TYPE(pw_r3d_rs_type)                                      :: virial_pw, v_laplace, v_laplacea, v_laplaceb
    1692              :       TYPE(section_vals_type), POINTER                   :: xc_fun_section
    1693              :       TYPE(xc_derivative_set_type)                       :: deriv_set1
    1694              :       TYPE(xc_rho_cflags_type)                           :: needs
    1695              :       TYPE(xc_rho_set_type)                              :: rho1_set, rho2_set
    1696              : 
    1697          340 :       CALL timeset(routineN, handle)
    1698              : 
    1699          340 :       my_calc_virial = .FALSE.
    1700          340 :       IF (PRESENT(calc_virial) .AND. PRESENT(virial_xc)) my_calc_virial = calc_virial
    1701              : 
    1702          340 :       nspins = SIZE(v_xc)
    1703              : 
    1704          340 :       NULLIFY (tau, tau_r, tau_a, tau_b)
    1705              : 
    1706          340 :       h = section_get_rval(xc_section, "STEP_SIZE")
    1707          340 :       nsteps = section_get_ival(xc_section, "NSTEPS")
    1708          340 :       IF (nsteps < LBOUND(weights, 2) .OR. nspins > UBOUND(weights, 2)) THEN
    1709            0 :          CPABORT("The number of steps must be a value from 1 to 4.")
    1710              :       END IF
    1711              : 
    1712          340 :       IF (nspins == 2) THEN
    1713          148 :          NULLIFY (vxc_rho, rho_g, vxc_tau)
    1714          444 :          ALLOCATE (rho_r(2))
    1715          444 :          DO ispin = 1, nspins
    1716          444 :             CALL pw_pool%create_pw(rho_r(ispin))
    1717              :          END DO
    1718          148 :          IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
    1719          162 :             ALLOCATE (tau_r(2))
    1720          162 :             DO ispin = 1, nspins
    1721          162 :                CALL pw_pool%create_pw(tau_r(ispin))
    1722              :             END DO
    1723              :          END IF
    1724          148 :          CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
    1725         1088 :          DO istep = -nsteps, nsteps
    1726          940 :             IF (istep == 0) CYCLE
    1727          792 :             weight = weights(istep, nsteps)/h
    1728          792 :             step = REAL(istep, dp)*h
    1729              :             CALL calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
    1730          792 :                                               tau_r, tau1_r, tau_a, tau_b, vxc_tau, xc_section, pw_pool, step)
    1731         2376 :             DO ispin = 1, nspins
    1732         1584 :                CALL pw_axpy(vxc_rho(ispin), v_xc(ispin), weight)
    1733         2376 :                IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
    1734          456 :                   CALL pw_axpy(vxc_tau(ispin), v_tau(ispin), weight)
    1735              :                END IF
    1736              :             END DO
    1737         2376 :             DO ispin = 1, nspins
    1738         2376 :                CALL vxc_rho(ispin)%release()
    1739              :             END DO
    1740          792 :             DEALLOCATE (vxc_rho)
    1741          940 :             IF (ASSOCIATED(vxc_tau)) THEN
    1742          684 :                DO ispin = 1, nspins
    1743          684 :                   CALL vxc_tau(ispin)%release()
    1744              :                END DO
    1745          228 :                DEALLOCATE (vxc_tau)
    1746              :             END IF
    1747              :          END DO
    1748          192 :       ELSE IF (nspins == 1 .AND. do_triplet) THEN
    1749           20 :          NULLIFY (vxc_rho, vxc_tau, rho_g)
    1750           60 :          ALLOCATE (rho_r(2))
    1751           60 :          DO ispin = 1, 2
    1752           60 :             CALL pw_pool%create_pw(rho_r(ispin))
    1753              :          END DO
    1754           20 :          IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
    1755            0 :             ALLOCATE (tau_r(2))
    1756            0 :             DO ispin = 1, nspins
    1757            0 :                CALL pw_pool%create_pw(tau_r(ispin))
    1758              :             END DO
    1759              :          END IF
    1760           20 :          CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rhoa=rhoa, rhob=rhob, tau_a=tau_a, tau_b=tau_b)
    1761          160 :          DO istep = -nsteps, nsteps
    1762          140 :             IF (istep == 0) CYCLE
    1763          120 :             weight = weights(istep, nsteps)/h
    1764          120 :             step = REAL(istep, dp)*h
    1765              :             ! K(alpha,alpha)
    1766          120 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
    1767              : !$OMP WORKSHARE
    1768              :             rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
    1769              : !$OMP END WORKSHARE NOWAIT
    1770              : !$OMP WORKSHARE
    1771              :             rho_r(2)%array(:, :, :) = rhob(:, :, :)
    1772              : !$OMP END WORKSHARE NOWAIT
    1773              :             IF (ASSOCIATED(tau1_r)) THEN
    1774              : !$OMP WORKSHARE
    1775              :                tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
    1776              : !$OMP END WORKSHARE NOWAIT
    1777              : !$OMP WORKSHARE
    1778              :                tau_r(2)%array(:, :, :) = tau_b(:, :, :)
    1779              : !$OMP END WORKSHARE NOWAIT
    1780              :             END IF
    1781              : !$OMP END PARALLEL
    1782              :             CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
    1783          120 :                                   pw_pool, .FALSE., virial_dummy)
    1784          120 :             CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
    1785          120 :             IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
    1786            0 :                CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
    1787              :             END IF
    1788          360 :             DO ispin = 1, 2
    1789          360 :                CALL vxc_rho(ispin)%release()
    1790              :             END DO
    1791          120 :             DEALLOCATE (vxc_rho)
    1792          120 :             IF (ASSOCIATED(vxc_tau)) THEN
    1793            0 :             DO ispin = 1, 2
    1794            0 :                CALL vxc_tau(ispin)%release()
    1795              :             END DO
    1796            0 :             DEALLOCATE (vxc_tau)
    1797              :             END IF
    1798          120 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
    1799              : !$OMP WORKSHARE
    1800              :             ! K(alpha,beta)
    1801              :             rho_r(1)%array(:, :, :) = rhoa(:, :, :)
    1802              : !$OMP END WORKSHARE NOWAIT
    1803              : !$OMP WORKSHARE
    1804              :             rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(1)%array(:, :, :)
    1805              : !$OMP END WORKSHARE NOWAIT
    1806              :             IF (ASSOCIATED(tau1_r)) THEN
    1807              : !$OMP WORKSHARE
    1808              :                tau_r(1)%array(:, :, :) = tau_a(:, :, :)
    1809              : !$OMP END WORKSHARE NOWAIT
    1810              : !$OMP WORKSHARE
    1811              :                tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(1)%array(:, :, :)
    1812              : !$OMP END WORKSHARE NOWAIT
    1813              :             END IF
    1814              : !$OMP END PARALLEL
    1815              :             CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
    1816          120 :                                   pw_pool, .FALSE., virial_dummy)
    1817          120 :             CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
    1818          120 :             IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
    1819            0 :                CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
    1820              :             END IF
    1821          360 :             DO ispin = 1, 2
    1822          360 :                CALL vxc_rho(ispin)%release()
    1823              :             END DO
    1824          120 :             DEALLOCATE (vxc_rho)
    1825          140 :             IF (ASSOCIATED(vxc_tau)) THEN
    1826            0 :             DO ispin = 1, 2
    1827            0 :                CALL vxc_tau(ispin)%release()
    1828              :             END DO
    1829            0 :             DEALLOCATE (vxc_tau)
    1830              :             END IF
    1831              :          END DO
    1832              :       ELSE
    1833          172 :          NULLIFY (vxc_rho, rho_r, rho_g, vxc_tau, tau_r, tau)
    1834          344 :          ALLOCATE (rho_r(1))
    1835          172 :          CALL pw_pool%create_pw(rho_r(1))
    1836          172 :          IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(v_tau)) THEN
    1837           92 :             ALLOCATE (tau_r(1))
    1838           46 :             CALL pw_pool%create_pw(tau_r(1))
    1839              :          END IF
    1840          172 :          CALL xc_rho_set_get(rho_set, can_return_null=.TRUE., rho=rho, tau=tau)
    1841         1280 :          DO istep = -nsteps, nsteps
    1842         1108 :             IF (istep == 0) CYCLE
    1843          936 :             weight = weights(istep, nsteps)/h
    1844          936 :             step = REAL(istep, dp)*h
    1845          936 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rho,step,rho1_r,tau1_r,tau,tau_r)
    1846              : !$OMP WORKSHARE
    1847              :             rho_r(1)%array(:, :, :) = rho(:, :, :) + step*rho1_r(1)%array(:, :, :)
    1848              : !$OMP END WORKSHARE NOWAIT
    1849              :             IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau) .AND. ASSOCIATED(tau1_r)) THEN
    1850              : !$OMP WORKSHARE
    1851              :                tau_r(1)%array(:, :, :) = tau(:, :, :) + step*tau1_r(1)%array(:, :, :)
    1852              : !$OMP END WORKSHARE NOWAIT
    1853              :             END IF
    1854              : !$OMP END PARALLEL
    1855              :             CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
    1856          936 :                                   pw_pool, .FALSE., virial_dummy)
    1857          936 :             CALL pw_axpy(vxc_rho(1), v_xc(1), weight)
    1858          936 :             IF (ASSOCIATED(vxc_tau) .AND. ASSOCIATED(v_tau)) THEN
    1859          276 :                CALL pw_axpy(vxc_tau(1), v_tau(1), weight)
    1860              :             END IF
    1861          936 :             CALL vxc_rho(1)%release()
    1862          936 :             DEALLOCATE (vxc_rho)
    1863         1108 :             IF (ASSOCIATED(vxc_tau)) THEN
    1864          276 :                CALL vxc_tau(1)%release()
    1865          276 :                DEALLOCATE (vxc_tau)
    1866              :             END IF
    1867              :          END DO
    1868              :       END IF
    1869              : 
    1870          340 :       IF (my_calc_virial) THEN
    1871           36 :          lsd = (nspins == 2)
    1872           36 :          IF (nspins == 1 .AND. do_triplet) THEN
    1873            0 :             lsd = .TRUE.
    1874              :          END IF
    1875              : 
    1876           36 :          CALL check_for_derivatives(deriv_set, (nspins == 2), rho_f, gradient_f, tau_f, laplace_f)
    1877              : 
    1878              :          ! Calculate the virial terms
    1879              :          ! Those arising from the first derivatives are treated like in xc_calc_2nd_deriv_analytical
    1880              :          ! Those arising from the second derivatives are calculated numerically
    1881              :          ! We assume that all metaGGA functionals require the gradient
    1882           36 :          IF (gradient_f) THEN
    1883          360 :             bo = rho_set%local_bounds
    1884              : 
    1885              :             ! Create the work grid for the virial terms
    1886           36 :             CALL allocate_pw(virial_pw, pw_pool, bo)
    1887              : 
    1888           36 :             gradient_cut = section_get_rval(xc_section, "GRADIENT_CUTOFF")
    1889              : 
    1890              :             ! create the container to store the argument of the functionals
    1891              :             CALL xc_rho_set_create(rho1_set, bo, &
    1892              :                                    rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
    1893              :                                    drho_cutoff=gradient_cut, &
    1894           36 :                                    tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
    1895              : 
    1896           36 :             xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
    1897           36 :             needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
    1898              : 
    1899              :             ! calculate the arguments needed by the functionals
    1900              :             CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
    1901              :                                    section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
    1902              :                                    section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
    1903           36 :                                    pw_pool)
    1904              : 
    1905           36 :             IF (lsd) THEN
    1906              :                CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, norm_drho=norm_drho, &
    1907              :                                    norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, &
    1908           10 :                                    laplace_rhoa=laplacea, laplace_rhob=laplaceb, can_return_null=.TRUE.)
    1909              :                CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, drhoa=drho1a, drhob=drho1b, laplace_rhoa=laplace1a, &
    1910           10 :                                    laplace_rhob=laplace1b, can_return_null=.TRUE.)
    1911              : 
    1912           10 :                CALL calc_drho_from_ab(drho, drhoa, drhob)
    1913           10 :                CALL calc_drho_from_ab(drho1, drho1a, drho1b)
    1914              :             ELSE
    1915           26 :                CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho, tau=tau, laplace_rho=laplace, can_return_null=.TRUE.)
    1916           26 :                CALL xc_rho_set_get(rho1_set, rho=rho1, drho=drho1, laplace_rho=laplace1, can_return_null=.TRUE.)
    1917              :             END IF
    1918              : 
    1919           36 :             CALL prepare_dr1dr(dr1dr, drho, drho1)
    1920              : 
    1921           36 :             IF (lsd) THEN
    1922           10 :                CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
    1923           10 :                CALL prepare_dr1dr(drb1drb, drhob, drho1b)
    1924              : 
    1925           10 :                CALL allocate_pw(v_drho, pw_pool, bo)
    1926           10 :                CALL allocate_pw(v_drhoa, pw_pool, bo)
    1927           10 :                CALL allocate_pw(v_drhob, pw_pool, bo)
    1928              : 
    1929           10 :                IF (ASSOCIATED(norm_drhoa)) CALL apply_drho(deriv_set, [deriv_norm_drhoa], virial_pw, drhoa, drho1a, virial_xc, &
    1930           10 :                                                            norm_drhoa, gradient_cut, dra1dra, v_drhoa%array)
    1931           10 :                IF (ASSOCIATED(norm_drhob)) CALL apply_drho(deriv_set, [deriv_norm_drhob], virial_pw, drhob, drho1b, virial_xc, &
    1932           10 :                                                            norm_drhob, gradient_cut, drb1drb, v_drhob%array)
    1933           10 :                IF (ASSOCIATED(norm_drho)) CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
    1934            6 :                                                           norm_drho, gradient_cut, dr1dr, v_drho%array)
    1935           10 :                IF (laplace_f) THEN
    1936            2 :                   CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa]), deriv_data=deriv_data)
    1937            2 :                   CPASSERT(ASSOCIATED(deriv_data))
    1938        15026 :                   virial_pw%array(:, :, :) = -rho1a(:, :, :)
    1939            2 :                   CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
    1940              : 
    1941            2 :                   CALL allocate_pw(v_laplacea, pw_pool, bo)
    1942              : 
    1943            2 :                   CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob]), deriv_data=deriv_data)
    1944            2 :                   CPASSERT(ASSOCIATED(deriv_data))
    1945        15026 :                   virial_pw%array(:, :, :) = -rho1b(:, :, :)
    1946            2 :                   CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
    1947              : 
    1948            2 :                   CALL allocate_pw(v_laplaceb, pw_pool, bo)
    1949              :                END IF
    1950              : 
    1951              :             ELSE
    1952              : 
    1953              :                ! Create the work grid for the potential of the gradient part
    1954           26 :                CALL allocate_pw(v_drho, pw_pool, bo)
    1955              : 
    1956              :                CALL apply_drho(deriv_set, [deriv_norm_drho], virial_pw, drho, drho1, virial_xc, &
    1957           26 :                                norm_drho, gradient_cut, dr1dr, v_drho%array)
    1958           26 :                IF (laplace_f) THEN
    1959            2 :                   CALL xc_derivative_get(xc_dset_get_derivative(deriv_set, [deriv_laplace_rho]), deriv_data=deriv_data)
    1960            2 :                   CPASSERT(ASSOCIATED(deriv_data))
    1961        28862 :                   virial_pw%array(:, :, :) = -rho1(:, :, :)
    1962            2 :                   CALL virial_laplace(virial_pw, pw_pool, virial_xc, deriv_data)
    1963              : 
    1964            2 :                   CALL allocate_pw(v_laplace, pw_pool, bo)
    1965              :                END IF
    1966              : 
    1967              :             END IF
    1968              : 
    1969           36 :             IF (lsd) THEN
    1970       150260 :                rho_r(1)%array = rhoa
    1971       150260 :                rho_r(2)%array = rhob
    1972              :             ELSE
    1973       701552 :                rho_r(1)%array = rho
    1974              :             END IF
    1975           36 :             IF (ASSOCIATED(tau1_r)) THEN
    1976            8 :             IF (lsd) THEN
    1977        60104 :                tau_r(1)%array = tau_a
    1978        60104 :                tau_r(2)%array = tau_b
    1979              :             ELSE
    1980       115448 :                tau_r(1)%array = tau
    1981              :             END IF
    1982              :             END IF
    1983              : 
    1984              :             ! Create deriv sets with same densities but different gradients
    1985           36 :             CALL xc_dset_create(deriv_set1, pw_pool)
    1986              : 
    1987           36 :             rho_cutoff = section_get_rval(xc_section, "DENSITY_CUTOFF")
    1988              : 
    1989              :             ! create the place where to store the argument for the functionals
    1990              :             CALL xc_rho_set_create(rho2_set, bo, &
    1991              :                                    rho_cutoff=rho_cutoff, &
    1992              :                                    drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
    1993           36 :                                    tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
    1994              : 
    1995              :             ! calculate the arguments needed by the functionals
    1996              :             CALL xc_rho_set_update(rho2_set, rho_r, rho_g, tau_r, needs, &
    1997              :                                    section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
    1998              :                                    section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
    1999           36 :                                    pw_pool)
    2000              : 
    2001           36 :             IF (lsd) THEN
    2002              :                CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b, tau_a=tau1a, tau_b=tau1b, &
    2003           10 :                                    laplace_rhoa=laplace1a, laplace_rhob=laplace1b, can_return_null=.TRUE.)
    2004              :                CALL xc_rho_set_get(rho2_set, norm_drhoa=norm_drho2a, norm_drhob=norm_drho2b, &
    2005           10 :                                    norm_drho=norm_drho2, laplace_rhoa=laplace2a, laplace_rhob=laplace2b, can_return_null=.TRUE.)
    2006              : 
    2007           64 :                DO istep = -nsteps, nsteps
    2008           54 :                   IF (istep == 0) CYCLE
    2009           44 :                   weight = weights(istep, nsteps)/h
    2010           44 :                   step = REAL(istep, dp)*h
    2011           44 :                   IF (ASSOCIATED(norm_drhoa)) THEN
    2012           44 :                      CALL get_derivs_rho(norm_drho2a, norm_drhoa, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2013              :                      CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
    2014           44 :                                            norm_drhoa, gradient_cut, weight, rho1a, v_drhoa%array)
    2015              :                      CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
    2016           44 :                                            norm_drhoa, gradient_cut, weight, rho1b, v_drhoa%array)
    2017              :                      CALL update_deriv_rho(deriv_set1, [deriv_norm_drhoa], bo, &
    2018           44 :                                            norm_drhoa, gradient_cut, weight, dra1dra, v_drhoa%array)
    2019              :                      CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
    2020           44 :                                                norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa%array, v_drhob%array)
    2021              :                      CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
    2022           44 :                                                norm_drhoa, gradient_cut, weight, dra1dra, dr1dr, v_drhoa%array, v_drho%array)
    2023           44 :                      IF (tau_f) THEN
    2024              :                         CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
    2025            8 :                                               norm_drhoa, gradient_cut, weight, tau1a, v_drhoa%array)
    2026              :                         CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
    2027            8 :                                               norm_drhoa, gradient_cut, weight, tau1b, v_drhoa%array)
    2028              :                      END IF
    2029           44 :                      IF (laplace_f) THEN
    2030              :                         CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
    2031            4 :                                               norm_drhoa, gradient_cut, weight, laplace1a, v_drhoa%array)
    2032              :                         CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
    2033            4 :                                               norm_drhoa, gradient_cut, weight, laplace1b, v_drhoa%array)
    2034              :                      END IF
    2035              :                   END IF
    2036              : 
    2037           44 :                   IF (ASSOCIATED(norm_drhob)) THEN
    2038           44 :                      CALL get_derivs_rho(norm_drho2b, norm_drhob, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2039              :                      CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
    2040           44 :                                            norm_drhob, gradient_cut, weight, rho1a, v_drhob%array)
    2041              :                      CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
    2042           44 :                                            norm_drhob, gradient_cut, weight, rho1b, v_drhob%array)
    2043              :                      CALL update_deriv_rho(deriv_set1, [deriv_norm_drhob], bo, &
    2044           44 :                                            norm_drhob, gradient_cut, weight, drb1drb, v_drhob%array)
    2045              :                      CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
    2046           44 :                                                norm_drhob, gradient_cut, weight, drb1drb, dra1dra, v_drhob%array, v_drhoa%array)
    2047              :                      CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drho], bo, &
    2048           44 :                                                norm_drhob, gradient_cut, weight, drb1drb, dr1dr, v_drhob%array, v_drho%array)
    2049           44 :                      IF (tau_f) THEN
    2050              :                         CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
    2051            8 :                                               norm_drhob, gradient_cut, weight, tau1a, v_drhob%array)
    2052              :                         CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
    2053            8 :                                               norm_drhob, gradient_cut, weight, tau1b, v_drhob%array)
    2054              :                      END IF
    2055           44 :                      IF (laplace_f) THEN
    2056              :                         CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
    2057            4 :                                               norm_drhob, gradient_cut, weight, laplace1a, v_drhob%array)
    2058              :                         CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
    2059            4 :                                               norm_drhob, gradient_cut, weight, laplace1b, v_drhob%array)
    2060              :                      END IF
    2061              :                   END IF
    2062              : 
    2063           44 :                   IF (ASSOCIATED(norm_drho)) THEN
    2064           20 :                      CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2065              :                      CALL update_deriv_rho(deriv_set1, [deriv_rhoa], bo, &
    2066           20 :                                            norm_drho, gradient_cut, weight, rho1a, v_drho%array)
    2067              :                      CALL update_deriv_rho(deriv_set1, [deriv_rhob], bo, &
    2068           20 :                                            norm_drho, gradient_cut, weight, rho1b, v_drho%array)
    2069              :                      CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
    2070           20 :                                            norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
    2071              :                      CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhoa], bo, &
    2072           20 :                                                norm_drho, gradient_cut, weight, dr1dr, dra1dra, v_drho%array, v_drhoa%array)
    2073              :                      CALL update_deriv_drho_ab(deriv_set1, [deriv_norm_drhob], bo, &
    2074           20 :                                                norm_drho, gradient_cut, weight, dr1dr, drb1drb, v_drho%array, v_drhob%array)
    2075           20 :                      IF (tau_f) THEN
    2076              :                         CALL update_deriv_rho(deriv_set1, [deriv_tau_a], bo, &
    2077            8 :                                               norm_drho, gradient_cut, weight, tau1a, v_drho%array)
    2078              :                         CALL update_deriv_rho(deriv_set1, [deriv_tau_b], bo, &
    2079            8 :                                               norm_drho, gradient_cut, weight, tau1b, v_drho%array)
    2080              :                      END IF
    2081           20 :                      IF (laplace_f) THEN
    2082              :                         CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhoa], bo, &
    2083            4 :                                               norm_drho, gradient_cut, weight, laplace1a, v_drho%array)
    2084              :                         CALL update_deriv_rho(deriv_set1, [deriv_laplace_rhob], bo, &
    2085            4 :                                               norm_drho, gradient_cut, weight, laplace1b, v_drho%array)
    2086              :                      END IF
    2087              :                   END IF
    2088              : 
    2089           54 :                   IF (laplace_f) THEN
    2090              : 
    2091            4 :                      CALL get_derivs_rho(laplace2a, laplacea, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2092              : 
    2093              :                      ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
    2094              :                      CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhoa], bo, &
    2095            4 :                                        weight, rho1a, v_laplacea%array)
    2096              :                      CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_rhob], bo, &
    2097            4 :                                        weight, rho1b, v_laplacea%array)
    2098            4 :                      IF (ASSOCIATED(norm_drho)) THEN
    2099              :                         CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drho], bo, &
    2100            4 :                                           weight, dr1dr, v_laplacea%array)
    2101              :                      END IF
    2102            4 :                      IF (ASSOCIATED(norm_drhoa)) THEN
    2103              :                         CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhoa], bo, &
    2104            4 :                                           weight, dra1dra, v_laplacea%array)
    2105              :                      END IF
    2106            4 :                      IF (ASSOCIATED(norm_drhob)) THEN
    2107              :                         CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_norm_drhob], bo, &
    2108            4 :                                           weight, drb1drb, v_laplacea%array)
    2109              :                      END IF
    2110              : 
    2111            4 :                      IF (ASSOCIATED(tau1a)) THEN
    2112              :                         CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_a], bo, &
    2113            4 :                                           weight, tau1a, v_laplacea%array)
    2114              :                      END IF
    2115            4 :                      IF (ASSOCIATED(tau1b)) THEN
    2116              :                         CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_tau_b], bo, &
    2117            4 :                                           weight, tau1b, v_laplacea%array)
    2118              :                      END IF
    2119              : 
    2120              :                      CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhoa], bo, &
    2121            4 :                                        weight, laplace1a, v_laplacea%array)
    2122              : 
    2123              :                      CALL update_deriv(deriv_set1, laplacea, rho_cutoff, [deriv_laplace_rhob], bo, &
    2124            4 :                                        weight, laplace1b, v_laplacea%array)
    2125              : 
    2126              :                      ! The same for the beta spin
    2127            4 :                      CALL get_derivs_rho(laplace2b, laplaceb, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2128              : 
    2129              :                      ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
    2130              :                      CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhoa], bo, &
    2131            4 :                                        weight, rho1a, v_laplaceb%array)
    2132              :                      CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_rhob], bo, &
    2133            4 :                                        weight, rho1b, v_laplaceb%array)
    2134            4 :                      IF (ASSOCIATED(norm_drho)) THEN
    2135              :                         CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drho], bo, &
    2136            4 :                                           weight, dr1dr, v_laplaceb%array)
    2137              :                      END IF
    2138            4 :                      IF (ASSOCIATED(norm_drhoa)) THEN
    2139              :                         CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhoa], bo, &
    2140            4 :                                           weight, dra1dra, v_laplaceb%array)
    2141              :                      END IF
    2142            4 :                      IF (ASSOCIATED(norm_drhob)) THEN
    2143              :                         CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_norm_drhob], bo, &
    2144            4 :                                           weight, drb1drb, v_laplaceb%array)
    2145              :                      END IF
    2146              : 
    2147            4 :                      IF (tau_f) THEN
    2148              :                         CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_a], bo, &
    2149            4 :                                           weight, tau1a, v_laplaceb%array)
    2150              :                         CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_tau_b], bo, &
    2151            4 :                                           weight, tau1b, v_laplaceb%array)
    2152              :                      END IF
    2153              : 
    2154              :                      CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhoa], bo, &
    2155            4 :                                        weight, laplace1a, v_laplaceb%array)
    2156              : 
    2157              :                      CALL update_deriv(deriv_set1, laplaceb, rho_cutoff, [deriv_laplace_rhob], bo, &
    2158            4 :                                        weight, laplace1b, v_laplaceb%array)
    2159              :                   END IF
    2160              :                END DO
    2161              : 
    2162           10 :                CALL virial_drho_drho(virial_pw, drhoa, v_drhoa, virial_xc)
    2163           10 :                CALL virial_drho_drho(virial_pw, drhob, v_drhob, virial_xc)
    2164           10 :                CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
    2165              : 
    2166           10 :                CALL deallocate_pw(v_drho, pw_pool)
    2167           10 :                CALL deallocate_pw(v_drhoa, pw_pool)
    2168           10 :                CALL deallocate_pw(v_drhob, pw_pool)
    2169              : 
    2170           10 :                IF (laplace_f) THEN
    2171        15026 :                   virial_pw%array(:, :, :) = -rhoa(:, :, :)
    2172            2 :                   CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplacea%array)
    2173            2 :                   CALL deallocate_pw(v_laplacea, pw_pool)
    2174              : 
    2175        15026 :                   virial_pw%array(:, :, :) = -rhob(:, :, :)
    2176            2 :                   CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplaceb%array)
    2177            2 :                   CALL deallocate_pw(v_laplaceb, pw_pool)
    2178              :                END IF
    2179              : 
    2180           10 :                CALL deallocate_pw(virial_pw, pw_pool)
    2181              : 
    2182           40 :                DO idir = 1, 3
    2183           30 :                   DEALLOCATE (drho(idir)%array)
    2184           40 :                   DEALLOCATE (drho1(idir)%array)
    2185              :                END DO
    2186           10 :                DEALLOCATE (dra1dra, drb1drb)
    2187              : 
    2188              :             ELSE
    2189           26 :                CALL xc_rho_set_get(rho1_set, rho=rho1, tau=tau1, laplace_rho=laplace1, can_return_null=.TRUE.)
    2190           26 :                CALL xc_rho_set_get(rho2_set, norm_drho=norm_drho2, laplace_rho=laplace2, can_return_null=.TRUE.)
    2191              : 
    2192          200 :                DO istep = -nsteps, nsteps
    2193          174 :                   IF (istep == 0) CYCLE
    2194          148 :                   weight = weights(istep, nsteps)/h
    2195          148 :                   step = REAL(istep, dp)*h
    2196          148 :                   CALL get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2197              : 
    2198              :                   ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
    2199              :                   CALL update_deriv_rho(deriv_set1, [deriv_rho], bo, &
    2200          148 :                                         norm_drho, gradient_cut, weight, rho1, v_drho%array)
    2201              :                   CALL update_deriv_rho(deriv_set1, [deriv_norm_drho], bo, &
    2202          148 :                                         norm_drho, gradient_cut, weight, dr1dr, v_drho%array)
    2203              : 
    2204          148 :                   IF (tau_f) THEN
    2205              :                      CALL update_deriv_rho(deriv_set1, [deriv_tau], bo, &
    2206           24 :                                            norm_drho, gradient_cut, weight, tau1, v_drho%array)
    2207              :                   END IF
    2208          174 :                   IF (laplace_f) THEN
    2209              :                      CALL update_deriv_rho(deriv_set1, [deriv_laplace_rho], bo, &
    2210           12 :                                            norm_drho, gradient_cut, weight, laplace1, v_drho%array)
    2211              : 
    2212           12 :                      CALL get_derivs_rho(laplace2, laplace, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2213              : 
    2214              :                      ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
    2215              :                      CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_rho], bo, &
    2216           12 :                                        weight, rho1, v_laplace%array)
    2217              :                      CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_norm_drho], bo, &
    2218           12 :                                        weight, dr1dr, v_laplace%array)
    2219              : 
    2220           12 :                      IF (tau_f) THEN
    2221              :                         CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_tau], bo, &
    2222           12 :                                           weight, tau1, v_laplace%array)
    2223              :                      END IF
    2224              : 
    2225              :                      CALL update_deriv(deriv_set1, laplace, rho_cutoff, [deriv_laplace_rho], bo, &
    2226           12 :                                        weight, laplace1, v_laplace%array)
    2227              :                   END IF
    2228              :                END DO
    2229              : 
    2230              :                ! Calculate the virial contribution from the potential
    2231           26 :                CALL virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
    2232              : 
    2233           26 :                CALL deallocate_pw(v_drho, pw_pool)
    2234              : 
    2235           26 :                IF (laplace_f) THEN
    2236        28862 :                   virial_pw%array(:, :, :) = -rho(:, :, :)
    2237            2 :                   CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace%array)
    2238            2 :                   CALL deallocate_pw(v_laplace, pw_pool)
    2239              :                END IF
    2240              : 
    2241           26 :                CALL deallocate_pw(virial_pw, pw_pool)
    2242              :             END IF
    2243              : 
    2244              :          END IF
    2245              : 
    2246           36 :          CALL xc_dset_release(deriv_set1)
    2247              : 
    2248           36 :          DEALLOCATE (dr1dr)
    2249              : 
    2250           36 :          CALL xc_rho_set_release(rho1_set)
    2251           36 :          CALL xc_rho_set_release(rho2_set)
    2252              :       END IF
    2253              : 
    2254          848 :       DO ispin = 1, SIZE(rho_r)
    2255          848 :          CALL pw_pool%give_back_pw(rho_r(ispin))
    2256              :       END DO
    2257          340 :       DEALLOCATE (rho_r)
    2258              : 
    2259          340 :       IF (ASSOCIATED(tau_r)) THEN
    2260          254 :       DO ispin = 1, SIZE(tau_r)
    2261          254 :          CALL pw_pool%give_back_pw(tau_r(ispin))
    2262              :       END DO
    2263          100 :       DEALLOCATE (tau_r)
    2264              :       END IF
    2265              : 
    2266          340 :       CALL timestop(handle)
    2267              : 
    2268        15300 :    END SUBROUTINE xc_calc_2nd_deriv_numerical
    2269              : 
    2270              : ! **************************************************************************************************
    2271              : !> \brief ...
    2272              : !> \param rho_r ...
    2273              : !> \param rho_g ...
    2274              : !> \param rho1_r ...
    2275              : !> \param rhoa ...
    2276              : !> \param rhob ...
    2277              : !> \param vxc_rho ...
    2278              : !> \param tau_r ...
    2279              : !> \param tau1_r ...
    2280              : !> \param tau_a ...
    2281              : !> \param tau_b ...
    2282              : !> \param vxc_tau ...
    2283              : !> \param xc_section ...
    2284              : !> \param pw_pool ...
    2285              : !> \param step ...
    2286              : ! **************************************************************************************************
    2287          792 :    SUBROUTINE calc_resp_potential_numer_ab(rho_r, rho_g, rho1_r, rhoa, rhob, vxc_rho, &
    2288              :                                            tau_r, tau1_r, tau_a, tau_b, vxc_tau, &
    2289              :                                            xc_section, pw_pool, step)
    2290              : 
    2291              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN) :: vxc_rho, vxc_tau
    2292              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN)            :: rho1_r
    2293              :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN), POINTER   :: tau1_r
    2294              :       TYPE(pw_pool_type), INTENT(IN), POINTER            :: pw_pool
    2295              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: xc_section
    2296              :       REAL(KIND=dp), INTENT(IN)                          :: step
    2297              :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER, INTENT(IN) :: rhoa, rhob, tau_a, tau_b
    2298              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN)   :: rho_r
    2299              :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
    2300              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               ::  tau_r
    2301              : 
    2302              :       CHARACTER(len=*), PARAMETER :: routineN = 'calc_resp_potential_numer_ab'
    2303              : 
    2304              :       INTEGER                                            :: handle
    2305              :       REAL(KIND=dp)                                      :: exc
    2306              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_dummy
    2307              : 
    2308          792 :       CALL timeset(routineN, handle)
    2309              : 
    2310          792 : !$OMP PARALLEL DEFAULT(NONE) SHARED(rho_r,rhoa,rhob,step,rho1_r,tau_r,tau_a,tau_b,tau1_r)
    2311              : !$OMP WORKSHARE
    2312              :       rho_r(1)%array(:, :, :) = rhoa(:, :, :) + step*rho1_r(1)%array(:, :, :)
    2313              : !$OMP END WORKSHARE NOWAIT
    2314              : !$OMP WORKSHARE
    2315              :       rho_r(2)%array(:, :, :) = rhob(:, :, :) + step*rho1_r(2)%array(:, :, :)
    2316              : !$OMP END WORKSHARE NOWAIT
    2317              :       IF (ASSOCIATED(tau1_r) .AND. ASSOCIATED(tau_r) .AND. ASSOCIATED(tau_a) .AND. ASSOCIATED(tau_b)) THEN
    2318              : !$OMP WORKSHARE
    2319              :          tau_r(1)%array(:, :, :) = tau_a(:, :, :) + step*tau1_r(1)%array(:, :, :)
    2320              : !$OMP END WORKSHARE NOWAIT
    2321              : !$OMP WORKSHARE
    2322              :          tau_r(2)%array(:, :, :) = tau_b(:, :, :) + step*tau1_r(2)%array(:, :, :)
    2323              : !$OMP END WORKSHARE NOWAIT
    2324              :       END IF
    2325              : !$OMP END PARALLEL
    2326              :       CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau_r, xc_section, &
    2327          792 :                             pw_pool, .FALSE., virial_dummy)
    2328              : 
    2329          792 :       CALL timestop(handle)
    2330              : 
    2331          792 :    END SUBROUTINE calc_resp_potential_numer_ab
    2332              : 
    2333              : ! **************************************************************************************************
    2334              : !> \brief calculates stress tensor and potential contributions from the first derivative
    2335              : !> \param deriv_set ...
    2336              : !> \param description ...
    2337              : !> \param virial_pw ...
    2338              : !> \param drho ...
    2339              : !> \param drho1 ...
    2340              : !> \param virial_xc ...
    2341              : !> \param norm_drho ...
    2342              : !> \param gradient_cut ...
    2343              : !> \param dr1dr ...
    2344              : !> \param v_drho ...
    2345              : ! **************************************************************************************************
    2346           52 :    SUBROUTINE apply_drho(deriv_set, description, virial_pw, drho, drho1, virial_xc, norm_drho, gradient_cut, dr1dr, v_drho)
    2347              : 
    2348              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
    2349              :       INTEGER, DIMENSION(:), INTENT(in)                  :: description
    2350              :       TYPE(pw_r3d_rs_type), INTENT(IN)                          :: virial_pw
    2351              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN)     :: drho, drho1
    2352              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: virial_xc
    2353              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: norm_drho
    2354              :       REAL(KIND=dp), INTENT(IN)                          :: gradient_cut
    2355              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: dr1dr
    2356              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: v_drho
    2357              : 
    2358              :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_drho'
    2359              : 
    2360              :       INTEGER                                            :: handle
    2361           52 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: deriv_data
    2362              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att
    2363              : 
    2364           52 :       CALL timeset(routineN, handle)
    2365              : 
    2366           52 :       deriv_att => xc_dset_get_derivative(deriv_set, description)
    2367           52 :       IF (ASSOCIATED(deriv_att)) THEN
    2368           52 :          CALL xc_derivative_get(deriv_att, deriv_data=deriv_data)
    2369           52 :          CALL virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
    2370              : 
    2371           52 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,gradient_cut,norm_drho,v_drho,deriv_data)
    2372              :          v_drho(:, :, :) = v_drho(:, :, :) + &
    2373              :                            deriv_data(:, :, :)*dr1dr(:, :, :)/MAX(gradient_cut, norm_drho(:, :, :))**2
    2374              : !$OMP END PARALLEL WORKSHARE
    2375              :       END IF
    2376              : 
    2377           52 :       CALL timestop(handle)
    2378              : 
    2379           52 :    END SUBROUTINE apply_drho
    2380              : 
    2381              : ! **************************************************************************************************
    2382              : !> \brief adds potential contributions from derivatives of rho or diagonal terms of norm_drho
    2383              : !> \param deriv_set1 ...
    2384              : !> \param description ...
    2385              : !> \param bo ...
    2386              : !> \param norm_drho norm_drho of which derivative is calculated
    2387              : !> \param gradient_cut ...
    2388              : !> \param h ...
    2389              : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
    2390              : !> \param v_drho ...
    2391              : ! **************************************************************************************************
    2392          728 :    SUBROUTINE update_deriv_rho(deriv_set1, description, bo, norm_drho, gradient_cut, weight, rho1, v_drho)
    2393              : 
    2394              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set1
    2395              :       INTEGER, DIMENSION(:), INTENT(in)                  :: description
    2396              :       INTEGER, DIMENSION(2, 3), INTENT(IN)               :: bo
    2397              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2398              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN)     :: norm_drho
    2399              :       REAL(KIND=dp), INTENT(IN)                          :: gradient_cut, weight
    2400              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2401              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN)     :: rho1
    2402              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2403              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT)  :: v_drho
    2404              : 
    2405              :       CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_rho'
    2406              : 
    2407              :       INTEGER                                            :: handle, i, j, k
    2408              :       REAL(KIND=dp)                                      :: de
    2409          728 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: deriv_data1
    2410              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att1
    2411              : 
    2412          728 :       CALL timeset(routineN, handle)
    2413              : 
    2414              :       ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
    2415          728 :       deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
    2416          728 :       IF (ASSOCIATED(deriv_att1)) THEN
    2417          728 :          CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
    2418              : !$OMP PARALLEL DO DEFAULT(NONE) &
    2419              : !$OMP             SHARED(bo,deriv_data1,weight,norm_drho,v_drho,rho1,gradient_cut) &
    2420              : !$OMP             PRIVATE(i,j,k,de) &
    2421          728 : !$OMP             COLLAPSE(3)
    2422              :          DO k = bo(1, 3), bo(2, 3)
    2423              :             DO j = bo(1, 2), bo(2, 2)
    2424              :                DO i = bo(1, 1), bo(2, 1)
    2425              :                   de = weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drho(i, j, k))**2
    2426              :                   v_drho(i, j, k) = v_drho(i, j, k) - de*rho1(i, j, k)
    2427              :                END DO
    2428              :             END DO
    2429              :          END DO
    2430              : !$OMP END PARALLEL DO
    2431              :       END IF
    2432              : 
    2433          728 :       CALL timestop(handle)
    2434              : 
    2435          728 :    END SUBROUTINE update_deriv_rho
    2436              : 
    2437              : ! **************************************************************************************************
    2438              : !> \brief adds potential contributions from derivatives of a component with positive and negative values
    2439              : !> \param deriv_set1 ...
    2440              : !> \param description ...
    2441              : !> \param bo ...
    2442              : !> \param h ...
    2443              : !> \param rho1 function to contract the derivative with (rho1 for rho, dr1dr for norm_drho)
    2444              : !> \param v ...
    2445              : ! **************************************************************************************************
    2446          120 :    SUBROUTINE update_deriv(deriv_set1, rho, rho_cutoff, description, bo, weight, rho1, v)
    2447              : 
    2448              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set1
    2449              :       INTEGER, DIMENSION(:), INTENT(in)                  :: description
    2450              :       INTEGER, DIMENSION(2, 3), INTENT(IN)               :: bo
    2451              :       REAL(KIND=dp), INTENT(IN)                          :: weight, rho_cutoff
    2452              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2453              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN)     :: rho, rho1
    2454              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2455              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT)  :: v
    2456              : 
    2457              :       CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv'
    2458              : 
    2459              :       INTEGER                                            :: handle, i, j, k
    2460              :       REAL(KIND=dp)                                      :: de
    2461          120 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: deriv_data1
    2462              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att1
    2463              : 
    2464          120 :       CALL timeset(routineN, handle)
    2465              : 
    2466              :       ! Obtain the numerical 2nd derivatives w.r.t. to drho and collect the potential
    2467          120 :       deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
    2468          120 :       IF (ASSOCIATED(deriv_att1)) THEN
    2469          120 :          CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
    2470              : !$OMP PARALLEL DO DEFAULT(NONE) &
    2471              : !$OMP             SHARED(bo,deriv_data1,weight,v,rho1,rho, rho_cutoff) &
    2472              : !$OMP             PRIVATE(i,j,k,de) &
    2473          120 : !$OMP             COLLAPSE(3)
    2474              :          DO k = bo(1, 3), bo(2, 3)
    2475              :             DO j = bo(1, 2), bo(2, 2)
    2476              :                DO i = bo(1, 1), bo(2, 1)
    2477              :                   ! We have to consider that the given density (mostly the Laplacian) may have positive and negative values
    2478              :                   de = weight*deriv_data1(i, j, k)/SIGN(MAX(ABS(rho(i, j, k)), rho_cutoff), rho(i, j, k))
    2479              :                   v(i, j, k) = v(i, j, k) + de*rho1(i, j, k)
    2480              :                END DO
    2481              :             END DO
    2482              :          END DO
    2483              : !$OMP END PARALLEL DO
    2484              :       END IF
    2485              : 
    2486          120 :       CALL timestop(handle)
    2487              : 
    2488          120 :    END SUBROUTINE update_deriv
    2489              : 
    2490              : ! **************************************************************************************************
    2491              : !> \brief adds mixed derivatives of norm_drho
    2492              : !> \param deriv_set1 ...
    2493              : !> \param description ...
    2494              : !> \param bo ...
    2495              : !> \param norm_drhoa norm_drho of which derivatives is calculated
    2496              : !> \param gradient_cut ...
    2497              : !> \param h ...
    2498              : !> \param dra1dra dr1dr corresponding to norm_drho
    2499              : !> \param drb1drb ...
    2500              : !> \param v_drhoa potential corresponding to norm_drho
    2501              : !> \param v_drhob ...
    2502              : ! **************************************************************************************************
    2503          216 :    SUBROUTINE update_deriv_drho_ab(deriv_set1, description, bo, &
    2504          216 :                                    norm_drhoa, gradient_cut, weight, dra1dra, drb1drb, v_drhoa, v_drhob)
    2505              : 
    2506              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set1
    2507              :       INTEGER, DIMENSION(:), INTENT(in)                  :: description
    2508              :       INTEGER, DIMENSION(2, 3), INTENT(IN)               :: bo
    2509              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2510              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN)     :: norm_drhoa
    2511              :       REAL(KIND=dp), INTENT(IN)                          :: gradient_cut, weight
    2512              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2513              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(IN)     :: dra1dra, drb1drb
    2514              :       REAL(KIND=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, &
    2515              :                                                      2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(INOUT)  :: v_drhoa, v_drhob
    2516              : 
    2517              :       CHARACTER(len=*), PARAMETER :: routineN = 'update_deriv_drho_ab'
    2518              : 
    2519              :       INTEGER                                            :: handle, i, j, k
    2520              :       REAL(KIND=dp)                                      :: de
    2521          216 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: deriv_data1
    2522              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att1
    2523              : 
    2524          216 :       CALL timeset(routineN, handle)
    2525              : 
    2526          216 :       deriv_att1 => xc_dset_get_derivative(deriv_set1, description)
    2527          216 :       IF (ASSOCIATED(deriv_att1)) THEN
    2528          168 :          CALL xc_derivative_get(deriv_att1, deriv_data=deriv_data1)
    2529              : !$OMP PARALLEL DO DEFAULT(NONE) &
    2530              : !$OMP             PRIVATE(k,j,i,de) &
    2531              : !$OMP             SHARED(bo,drb1drb,dra1dra,deriv_data1,weight,gradient_cut,norm_drhoa,v_drhoa,v_drhob) &
    2532          168 : !$OMP             COLLAPSE(3)
    2533              :          DO k = bo(1, 3), bo(2, 3)
    2534              :             DO j = bo(1, 2), bo(2, 2)
    2535              :                DO i = bo(1, 1), bo(2, 1)
    2536              :                   ! We introduce a factor of two because we will average between both numerical derivatives
    2537              :                   de = 0.5_dp*weight*deriv_data1(i, j, k)/MAX(gradient_cut, norm_drhoa(i, j, k))**2
    2538              :                   v_drhoa(i, j, k) = v_drhoa(i, j, k) - de*drb1drb(i, j, k)
    2539              :                   v_drhob(i, j, k) = v_drhob(i, j, k) - de*dra1dra(i, j, k)
    2540              :                END DO
    2541              :             END DO
    2542              :          END DO
    2543              : !$OMP END PARALLEL DO
    2544              :       END IF
    2545              : 
    2546          216 :       CALL timestop(handle)
    2547              : 
    2548          216 :    END SUBROUTINE update_deriv_drho_ab
    2549              : 
    2550              : ! **************************************************************************************************
    2551              : !> \brief calculate derivative sets for helper points
    2552              : !> \param norm_drho2 norm_drho of new points
    2553              : !> \param norm_drho norm_drho of KS density
    2554              : !> \param h ...
    2555              : !> \param xc_fun_section ...
    2556              : !> \param lsd ...
    2557              : !> \param rho2_set rho_set for new points
    2558              : !> \param deriv_set1 will contain derivatives of the perturbed density
    2559              : ! **************************************************************************************************
    2560          276 :    SUBROUTINE get_derivs_rho(norm_drho2, norm_drho, step, xc_fun_section, lsd, rho2_set, deriv_set1)
    2561              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT)     :: norm_drho2
    2562              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: norm_drho
    2563              :       REAL(KIND=dp), INTENT(IN)                          :: step
    2564              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: xc_fun_section
    2565              :       LOGICAL, INTENT(IN)                                :: lsd
    2566              :       TYPE(xc_rho_set_type), INTENT(INOUT)               :: rho2_set
    2567              :       TYPE(xc_derivative_set_type)                       :: deriv_set1
    2568              : 
    2569              :       CHARACTER(len=*), PARAMETER :: routineN = 'get_derivs_rho'
    2570              : 
    2571              :       INTEGER                                            :: handle
    2572              : 
    2573          276 :       CALL timeset(routineN, handle)
    2574              : 
    2575              :       ! Copy the densities, do one step into the direction of drho
    2576          276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2,step)
    2577              :       norm_drho2 = norm_drho*(1.0_dp + step)
    2578              : !$OMP END PARALLEL WORKSHARE
    2579              : 
    2580          276 :       CALL xc_dset_zero_all(deriv_set1)
    2581              : 
    2582              :       ! Calculate the derivatives of the functional
    2583              :       CALL xc_functionals_eval(xc_fun_section, &
    2584              :                                lsd=lsd, &
    2585              :                                rho_set=rho2_set, &
    2586              :                                deriv_set=deriv_set1, &
    2587          276 :                                deriv_order=1)
    2588              : 
    2589              :       ! Return to the original values
    2590          276 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(norm_drho,norm_drho2)
    2591              :       norm_drho2 = norm_drho
    2592              : !$OMP END PARALLEL WORKSHARE
    2593              : 
    2594          276 :       CALL divide_by_norm_drho(deriv_set1, rho2_set, lsd)
    2595              : 
    2596          276 :       CALL timestop(handle)
    2597              : 
    2598          276 :    END SUBROUTINE get_derivs_rho
    2599              : 
    2600              : ! **************************************************************************************************
    2601              : !> \brief Calculates the second derivative of E_xc at rho in the direction
    2602              : !>      rho1  (if you see the second derivative as bilinear form)
    2603              : !>      partial_rho|_(rho=rho) partial_rho|_(rho=rho) E_xc drho(rho1)drho
    2604              : !>      The other direction is still undetermined, thus it returns
    2605              : !>      a potential (partial integration is performed to reduce it to
    2606              : !>      function of rho, removing the dependence from its partial derivs)
    2607              : !>      Has to be called after the setup by xc_prep_2nd_deriv.
    2608              : !> \param v_xc       exchange-correlation potential
    2609              : !> \param v_xc_tau ...
    2610              : !> \param deriv_set  derivatives of the exchange-correlation potential
    2611              : !> \param rho_set    object containing the density at which the derivatives were calculated
    2612              : !> \param rho1_set   object containing the density with which to fold
    2613              : !> \param pw_pool    the pool for the grids
    2614              : !> \param xc_section XC parameters
    2615              : !> \param gapw       Gaussian and augmented plane waves calculation
    2616              : !> \param vxg ...
    2617              : !> \param tddfpt_fac factor that multiplies the crossterms (tddfpt triplets
    2618              : !>        on a closed shell system it should be -1, defaults to 1)
    2619              : !> \param compute_virial ...
    2620              : !> \param virial_xc ...
    2621              : !> \note
    2622              : !>      The old version of this routine was smarter: it handled split_desc(1)
    2623              : !>      and split_desc(2) separately, thus the code automatically handled all
    2624              : !>      possible cross terms (you only had to check if it was diagonal to avoid
    2625              : !>      double counting). I think that is the way to go if you want to add more
    2626              : !>      terms (tau,rho in LSD,...). The problem with the old code was that it
    2627              : !>      because of the old functional structure it sometime guessed wrongly
    2628              : !>      which derivative was where. There were probably still bugs with gradient
    2629              : !>      corrected functionals (never tested), and it didn't contain first
    2630              : !>      derivatives with respect to drho (that contribute also to the second
    2631              : !>      derivative wrt. rho).
    2632              : !>      The code was a little complex because it really tried to handle any
    2633              : !>      functional derivative in the most efficient way with the given contents of
    2634              : !>      rho_set.
    2635              : !>      Anyway I strongly encourage whoever wants to modify this code to give a
    2636              : !>      look to the old version. [fawzi]
    2637              : ! **************************************************************************************************
    2638        30762 :    SUBROUTINE xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, &
    2639              :                                            pw_pool, xc_section, gapw, vxg, tddfpt_fac, &
    2640              :                                            compute_virial, virial_xc)
    2641              : 
    2642              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER             :: v_xc, v_xc_tau
    2643              :       TYPE(xc_derivative_set_type)                       :: deriv_set
    2644              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set, rho1_set
    2645              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
    2646              :       TYPE(section_vals_type), POINTER                   :: xc_section
    2647              :       LOGICAL, INTENT(IN), OPTIONAL                      :: gapw
    2648              :       REAL(kind=dp), DIMENSION(:, :, :, :), OPTIONAL, &
    2649              :          POINTER                                         :: vxg
    2650              :       REAL(kind=dp), INTENT(in), OPTIONAL                :: tddfpt_fac
    2651              :       LOGICAL, INTENT(IN), OPTIONAL                      :: compute_virial
    2652              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
    2653              :          OPTIONAL                                        :: virial_xc
    2654              : 
    2655              :       CHARACTER(len=*), PARAMETER :: routineN = 'xc_calc_2nd_deriv_analytical'
    2656              : 
    2657              :       INTEGER                                            :: handle, i, ia, idir, ir, ispin, j, jdir, &
    2658              :                                                             k, nspins, xc_deriv_method_id
    2659              :       INTEGER, DIMENSION(2, 3)                           :: bo
    2660              :       LOGICAL                                            :: gradient_f, lsd, my_compute_virial, &
    2661              :                                                             my_gapw, tau_f, laplace_f, rho_f
    2662              :       REAL(KIND=dp)                                      :: fac, gradient_cut, tmp, factor2
    2663        30762 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dr1dr, dra1dra, drb1drb
    2664        61524 :       REAL(kind=dp), DIMENSION(:, :, :), POINTER         :: deriv_data, e_drhoa, e_drhob, &
    2665        30762 :                                                             e_drho, norm_drho, norm_drhoa, &
    2666        30762 :                                                             norm_drhob, rho1, rho1a, rho1b, &
    2667        30762 :                                                             tau1, tau1a, tau1b, laplace1, laplace1a, laplace1b, &
    2668        30762 :                                                             rho, rhoa, rhob
    2669       584478 :       TYPE(cp_3d_r_cp_type), DIMENSION(3)                :: drho, drho1, drho1a, drho1b, drhoa, drhob
    2670        30762 :       TYPE(pw_r3d_rs_type), DIMENSION(:), ALLOCATABLE         :: v_drhoa, v_drhob, v_drho, v_laplace
    2671        30762 :       TYPE(pw_r3d_rs_type), DIMENSION(:, :), ALLOCATABLE      :: v_drho_r
    2672              :       TYPE(pw_r3d_rs_type)                                      ::  virial_pw
    2673              :       TYPE(pw_c1d_gs_type) :: tmp_g, vxc_g
    2674              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att
    2675              : 
    2676        30762 :       CALL timeset(routineN, handle)
    2677              : 
    2678        30762 :       NULLIFY (e_drhoa, e_drhob, e_drho)
    2679              : 
    2680        30762 :       my_gapw = .FALSE.
    2681        30762 :       IF (PRESENT(gapw)) my_gapw = gapw
    2682              : 
    2683        30762 :       my_compute_virial = .FALSE.
    2684        30762 :       IF (PRESENT(compute_virial)) my_compute_virial = compute_virial
    2685              : 
    2686        30762 :       CPASSERT(ASSOCIATED(v_xc))
    2687        30762 :       CPASSERT(ASSOCIATED(xc_section))
    2688        30762 :       IF (my_gapw) THEN
    2689         9842 :          CPASSERT(PRESENT(vxg))
    2690              :       END IF
    2691        30762 :       IF (my_compute_virial) THEN
    2692          358 :          CPASSERT(PRESENT(virial_xc))
    2693              :       END IF
    2694              : 
    2695              :       CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
    2696        30762 :                                 i_val=xc_deriv_method_id)
    2697        30762 :       CALL xc_rho_set_get(rho_set, drho_cutoff=gradient_cut)
    2698        30762 :       nspins = SIZE(v_xc)
    2699        30762 :       lsd = ASSOCIATED(rho_set%rhoa)
    2700        30762 :       fac = 0.0_dp
    2701        30762 :       factor2 = 1.0_dp
    2702        30762 :       IF (PRESENT(tddfpt_fac)) fac = tddfpt_fac
    2703        30762 :       IF (PRESENT(tddfpt_fac)) factor2 = tddfpt_fac
    2704              : 
    2705       307620 :       bo = rho_set%local_bounds
    2706              : 
    2707        30762 :       CALL check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
    2708              : 
    2709        30762 :       IF (tau_f) THEN
    2710          656 :          CPASSERT(ASSOCIATED(v_xc_tau))
    2711              :       END IF
    2712              : 
    2713        30762 :       IF (gradient_f) THEN
    2714       198850 :          ALLOCATE (v_drho_r(3, nspins), v_drho(nspins))
    2715        39770 :          DO ispin = 1, nspins
    2716        82120 :             DO idir = 1, 3
    2717        82120 :                CALL allocate_pw(v_drho_r(idir, ispin), pw_pool, bo)
    2718              :             END DO
    2719        39770 :             CALL allocate_pw(v_drho(ispin), pw_pool, bo)
    2720              :          END DO
    2721              : 
    2722        19240 :          IF (xc_requires_tmp_g(xc_deriv_method_id) .AND. .NOT. my_gapw) THEN
    2723        11576 :             IF (ASSOCIATED(pw_pool)) THEN
    2724        11576 :                CALL pw_pool%create_pw(tmp_g)
    2725        11576 :                CALL pw_pool%create_pw(vxc_g)
    2726              :             ELSE
    2727              :                ! remember to refix for gapw
    2728            0 :                CPABORT("XC_DERIV method is not implemented in GAPW")
    2729              :             END IF
    2730              :          END IF
    2731              :       END IF
    2732              : 
    2733        64718 :       DO ispin = 1, nspins
    2734    958202225 :          v_xc(ispin)%array = 0.0_dp
    2735              :       END DO
    2736              : 
    2737        30762 :       IF (tau_f) THEN
    2738         1386 :          DO ispin = 1, nspins
    2739     14897758 :             v_xc_tau(ispin)%array = 0.0_dp
    2740              :          END DO
    2741              :       END IF
    2742              : 
    2743        30762 :       IF (laplace_f .AND. my_gapw) &
    2744            0 :          CPABORT("Laplace-dependent functional not implemented with GAPW!")
    2745              : 
    2746        30762 :       IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) CALL allocate_pw(virial_pw, pw_pool, bo)
    2747              : 
    2748        30762 :       IF (lsd) THEN
    2749              : 
    2750              :          !-------------------!
    2751              :          ! UNrestricted case !
    2752              :          !-------------------!
    2753              : 
    2754         4232 :          CALL xc_rho_set_get(rho1_set, rhoa=rho1a, rhob=rho1b)
    2755              : 
    2756         4232 :          IF (gradient_f) THEN
    2757              :             CALL xc_rho_set_get(rho_set, drhoa=drhoa, drhob=drhob, &
    2758         2000 :                                 norm_drho=norm_drho, norm_drhoa=norm_drhoa, norm_drhob=norm_drhob)
    2759         2000 :             CALL xc_rho_set_get(rho1_set, drhoa=drho1a, drhob=drho1b)
    2760              : 
    2761         2000 :             CALL calc_drho_from_ab(drho, drhoa, drhob)
    2762         2000 :             CALL calc_drho_from_ab(drho1, drho1a, drho1b)
    2763              : 
    2764         2000 :             CALL prepare_dr1dr(dra1dra, drhoa, drho1a)
    2765         2000 :             IF (nspins /= 1) THEN
    2766         1290 :                CALL prepare_dr1dr(drb1drb, drhob, drho1b)
    2767         1290 :                CALL prepare_dr1dr(dr1dr, drho, drho1)
    2768              :             ELSE
    2769          710 :                CALL prepare_dr1dr(drb1drb, drhob, drho1b)
    2770          710 :                CALL prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
    2771              :             END IF
    2772              : 
    2773        14580 :             ALLOCATE (v_drhoa(nspins), v_drhob(nspins))
    2774         5290 :             DO ispin = 1, nspins
    2775         3290 :                CALL allocate_pw(v_drhoa(ispin), pw_pool, bo)
    2776         5290 :                CALL allocate_pw(v_drhob(ispin), pw_pool, bo)
    2777              :             END DO
    2778              : 
    2779              :          END IF
    2780              : 
    2781         4232 :          IF (laplace_f) THEN
    2782           38 :             CALL xc_rho_set_get(rho1_set, laplace_rhoa=laplace1a, laplace_rhob=laplace1b)
    2783              : 
    2784          190 :             ALLOCATE (v_laplace(nspins))
    2785          114 :             DO ispin = 1, nspins
    2786          114 :                CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
    2787              :             END DO
    2788              : 
    2789           38 :             IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob)
    2790              :          END IF
    2791              : 
    2792         4232 :          IF (tau_f) THEN
    2793           74 :             CALL xc_rho_set_get(rho1_set, tau_a=tau1a, tau_b=tau1b)
    2794              :          END IF
    2795              : 
    2796         4232 :          IF (nspins /= 1) THEN
    2797              : 
    2798        33242 :             $:add_2nd_derivative_terms(arguments_openshell)
    2799              : 
    2800              :          ELSE
    2801              : 
    2802         1038 :             $:add_2nd_derivative_terms(arguments_triplet_outer, arguments_triplet_inner)
    2803              : 
    2804              :          END IF
    2805              : 
    2806         4232 :          IF (gradient_f) THEN
    2807              : 
    2808         2000 :             IF (my_compute_virial) THEN
    2809           10 :                CALL virial_drho_drho(virial_pw, drhoa, v_drhoa(1), virial_xc)
    2810           10 :                CALL virial_drho_drho(virial_pw, drhob, v_drhob(2), virial_xc)
    2811           40 :                DO idir = 1, 3
    2812           30 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
    2813              :                   virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*(v_drho(1)%array(:, :, :) + v_drho(2)%array(:, :, :))
    2814              : !$OMP END PARALLEL WORKSHARE
    2815          100 :                   DO jdir = 1, idir
    2816              :                      tmp = -0.5_dp*virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
    2817           60 :                                                                                drho(jdir)%array(:, :, :))
    2818           60 :                      virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
    2819           90 :                      virial_xc(idir, jdir) = virial_xc(jdir, idir)
    2820              :                   END DO
    2821              :                END DO
    2822              :             END IF ! my_compute_virial
    2823              : 
    2824         2000 :             IF (my_gapw) THEN
    2825              : !$OMP PARALLEL DO DEFAULT(NONE) &
    2826              : !$OMP             PRIVATE(ia,idir,ispin,ir) &
    2827              : !$OMP             SHARED(bo,nspins,vxg,drhoa,drhob,v_drhoa,v_drhob,v_drho, &
    2828              : !$OMP                    e_drhoa,e_drhob,e_drho,drho1a,drho1b,fac,drho,drho1) &
    2829          524 : !$OMP             COLLAPSE(3)
    2830              :                DO ir = bo(1, 2), bo(2, 2)
    2831              :                   DO ia = bo(1, 1), bo(2, 1)
    2832              :                      DO idir = 1, 3
    2833              :                         DO ispin = 1, nspins
    2834              :                            vxg(idir, ia, ir, ispin) = &
    2835              :                               -(v_drhoa(ispin)%array(ia, ir, 1)*drhoa(idir)%array(ia, ir, 1) + &
    2836              :                                 v_drhob(ispin)%array(ia, ir, 1)*drhob(idir)%array(ia, ir, 1) + &
    2837              :                                 v_drho(ispin)%array(ia, ir, 1)*drho(idir)%array(ia, ir, 1))
    2838              :                         END DO
    2839              :                         IF (ASSOCIATED(e_drhoa)) THEN
    2840              :                            vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
    2841              :                                                   e_drhoa(ia, ir, 1)*drho1a(idir)%array(ia, ir, 1)
    2842              :                         END IF
    2843              :                         IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
    2844              :                            vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
    2845              :                                                   e_drhob(ia, ir, 1)*drho1b(idir)%array(ia, ir, 1)
    2846              :                         END IF
    2847              :                         IF (ASSOCIATED(e_drho)) THEN
    2848              :                            IF (nspins /= 1) THEN
    2849              :                               vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
    2850              :                                                      e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
    2851              :                               vxg(idir, ia, ir, 2) = vxg(idir, ia, ir, 2) + &
    2852              :                                                      e_drho(ia, ir, 1)*drho1(idir)%array(ia, ir, 1)
    2853              :                            ELSE
    2854              :                               vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + &
    2855              :                                                      e_drho(ia, ir, 1)*(drho1a(idir)%array(ia, ir, 1) + &
    2856              :                                                                         fac*drho1b(idir)%array(ia, ir, 1))
    2857              :                            END IF
    2858              :                         END IF
    2859              :                      END DO
    2860              :                   END DO
    2861              :                END DO
    2862              : !$OMP END PARALLEL DO
    2863              :             ELSE
    2864              : 
    2865              :                ! partial integration
    2866         5904 :                DO idir = 1, 3
    2867              : 
    2868        12342 :                   DO ispin = 1, nspins
    2869        12342 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,v_drhoa,v_drhob,v_drho,drhoa,drhob,drho,ispin,idir)
    2870              :                      v_drho_r(idir, ispin)%array(:, :, :) = &
    2871              :                         v_drhoa(ispin)%array(:, :, :)*drhoa(idir)%array(:, :, :) + &
    2872              :                         v_drhob(ispin)%array(:, :, :)*drhob(idir)%array(:, :, :) + &
    2873              :                         v_drho(ispin)%array(:, :, :)*drho(idir)%array(:, :, :)
    2874              : !$OMP END PARALLEL WORKSHARE
    2875              :                   END DO
    2876         4428 :                   IF (ASSOCIATED(e_drhoa)) THEN
    2877         4428 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,e_drhoa,drho1a,idir)
    2878              :                      v_drho_r(idir, 1)%array(:, :, :) = v_drho_r(idir, 1)%array(:, :, :) - &
    2879              :                                                         e_drhoa(:, :, :)*drho1a(idir)%array(:, :, :)
    2880              : !$OMP END PARALLEL WORKSHARE
    2881              :                   END IF
    2882         4428 :                   IF (nspins /= 1 .AND. ASSOCIATED(e_drhob)) THEN
    2883         3486 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,e_drhob,drho1b,idir)
    2884              :                      v_drho_r(idir, 2)%array(:, :, :) = v_drho_r(idir, 2)%array(:, :, :) - &
    2885              :                                                         e_drhob(:, :, :)*drho1b(idir)%array(:, :, :)
    2886              : !$OMP END PARALLEL WORKSHARE
    2887              :                   END IF
    2888         5904 :                   IF (ASSOCIATED(e_drho)) THEN
    2889              : !$OMP PARALLEL DO DEFAULT(NONE) &
    2890              : !$OMP              PRIVATE(k,j,i) &
    2891              : !$OMP              SHARED(bo,v_drho_r,e_drho,drho1a,drho1b,drho1,fac,idir,nspins) &
    2892         4428 : !$OMP              COLLAPSE(3)
    2893              :                      DO k = bo(1, 3), bo(2, 3)
    2894              :                         DO j = bo(1, 2), bo(2, 2)
    2895              :                            DO i = bo(1, 1), bo(2, 1)
    2896              :                               IF (nspins /= 1) THEN
    2897              :                                  v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
    2898              :                                                                     e_drho(i, j, k)*drho1(idir)%array(i, j, k)
    2899              :                                  v_drho_r(idir, 2)%array(i, j, k) = v_drho_r(idir, 2)%array(i, j, k) - &
    2900              :                                                                     e_drho(i, j, k)*drho1(idir)%array(i, j, k)
    2901              :                               ELSE
    2902              :                                  v_drho_r(idir, 1)%array(i, j, k) = v_drho_r(idir, 1)%array(i, j, k) - &
    2903              :                                                                     e_drho(i, j, k)*(drho1a(idir)%array(i, j, k) + &
    2904              :                                                                                      fac*drho1b(idir)%array(i, j, k))
    2905              :                               END IF
    2906              :                            END DO
    2907              :                         END DO
    2908              :                      END DO
    2909              : !$OMP END PARALLEL DO
    2910              :                   END IF
    2911              :                END DO
    2912              : 
    2913         4114 :                DO ispin = 1, nspins
    2914              :                   ! partial integration
    2915         4114 :                   CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, ispin), tmp_g, vxc_g, v_xc(ispin))
    2916              :                END DO ! ispin
    2917              : 
    2918              :             END IF
    2919              : 
    2920         8000 :             DO idir = 1, 3
    2921         6000 :                DEALLOCATE (drho(idir)%array)
    2922         8000 :                DEALLOCATE (drho1(idir)%array)
    2923              :             END DO
    2924              : 
    2925         5290 :             DO ispin = 1, nspins
    2926         3290 :                CALL deallocate_pw(v_drhoa(ispin), pw_pool)
    2927         5290 :                CALL deallocate_pw(v_drhob(ispin), pw_pool)
    2928              :             END DO
    2929              : 
    2930         2000 :             DEALLOCATE (v_drhoa, v_drhob)
    2931              : 
    2932              :          END IF ! gradient_f
    2933              : 
    2934         4232 :          IF (laplace_f .AND. my_compute_virial) THEN
    2935        15026 :             virial_pw%array(:, :, :) = -rhoa(:, :, :)
    2936            2 :             CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
    2937        15026 :             virial_pw%array(:, :, :) = -rhob(:, :, :)
    2938            2 :             CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(2)%array)
    2939              :          END IF
    2940              : 
    2941              :       ELSE
    2942              : 
    2943              :          !-----------------!
    2944              :          ! restricted case !
    2945              :          !-----------------!
    2946              : 
    2947        26530 :          CALL xc_rho_set_get(rho1_set, rho=rho1)
    2948              : 
    2949        26530 :          IF (gradient_f) THEN
    2950        17240 :             CALL xc_rho_set_get(rho_set, drho=drho, norm_drho=norm_drho)
    2951        17240 :             CALL xc_rho_set_get(rho1_set, drho=drho1)
    2952        17240 :             CALL prepare_dr1dr(dr1dr, drho, drho1)
    2953              :          END IF
    2954              : 
    2955        26530 :          IF (laplace_f) THEN
    2956          136 :             CALL xc_rho_set_get(rho1_set, laplace_rho=laplace1)
    2957              : 
    2958          544 :             ALLOCATE (v_laplace(nspins))
    2959          272 :             DO ispin = 1, nspins
    2960          272 :                CALL allocate_pw(v_laplace(ispin), pw_pool, bo)
    2961              :             END DO
    2962              : 
    2963          136 :             IF (my_compute_virial) CALL xc_rho_set_get(rho_set, rho=rho)
    2964              :          END IF
    2965              : 
    2966        26530 :          IF (tau_f) THEN
    2967          582 :             CALL xc_rho_set_get(rho1_set, tau=tau1)
    2968              :          END IF
    2969              : 
    2970       321046 :          $:add_2nd_derivative_terms(arguments_closedshell)
    2971              : 
    2972        26530 :          IF (gradient_f) THEN
    2973              : 
    2974        17240 :             IF (my_compute_virial) THEN
    2975          222 :                CALL virial_drho_drho(virial_pw, drho, v_drho(1), virial_xc)
    2976              :             END IF ! my_compute_virial
    2977              : 
    2978        17240 :             IF (my_gapw) THEN
    2979              : 
    2980        26896 :                DO idir = 1, 3
    2981              : !$OMP PARALLEL DO DEFAULT(NONE) &
    2982              : !$OMP             PRIVATE(ia,ir) &
    2983              : !$OMP             SHARED(bo,vxg,drho,v_drho,e_drho,drho1,idir,factor2) &
    2984        26896 : !$OMP             COLLAPSE(2)
    2985              :                   DO ia = bo(1, 1), bo(2, 1)
    2986              :                      DO ir = bo(1, 2), bo(2, 2)
    2987              :                         vxg(idir, ia, ir, 1) = -drho(idir)%array(ia, ir, 1)*v_drho(1)%array(ia, ir, 1)
    2988              :                         IF (ASSOCIATED(e_drho)) THEN
    2989              :                            vxg(idir, ia, ir, 1) = vxg(idir, ia, ir, 1) + factor2*drho1(idir)%array(ia, ir, 1)*e_drho(ia, ir, 1)
    2990              :                         END IF
    2991              :                      END DO
    2992              :                   END DO
    2993              : !$OMP END PARALLEL DO
    2994              :                END DO
    2995              : 
    2996              :             ELSE
    2997              :                ! partial integration
    2998        42064 :                DO idir = 1, 3
    2999        42064 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(v_drho_r,drho,v_drho,drho1,e_drho,idir)
    3000              :                   v_drho_r(idir, 1)%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho(1)%array(:, :, :) - &
    3001              :                                                      drho1(idir)%array(:, :, :)*e_drho(:, :, :)
    3002              : !$OMP END PARALLEL WORKSHARE
    3003              :                END DO
    3004              : 
    3005        10516 :                CALL xc_pw_divergence(xc_deriv_method_id, v_drho_r(:, 1), tmp_g, vxc_g, v_xc(1))
    3006              :             END IF
    3007              : 
    3008              :          END IF
    3009              : 
    3010        26530 :          IF (laplace_f .AND. my_compute_virial) THEN
    3011       294530 :             virial_pw%array(:, :, :) = -rho(:, :, :)
    3012           14 :             CALL virial_laplace(virial_pw, pw_pool, virial_xc, v_laplace(1)%array)
    3013              :          END IF
    3014              : 
    3015              :       END IF
    3016              : 
    3017        30762 :       IF (laplace_f) THEN
    3018          386 :          DO ispin = 1, nspins
    3019          212 :             CALL xc_pw_laplace(v_laplace(ispin), pw_pool, xc_deriv_method_id)
    3020          386 :             CALL pw_axpy(v_laplace(ispin), v_xc(ispin))
    3021              :          END DO
    3022              :       END IF
    3023              : 
    3024        30762 :       IF (gradient_f) THEN
    3025              : 
    3026        39770 :          DO ispin = 1, nspins
    3027        20530 :             CALL deallocate_pw(v_drho(ispin), pw_pool)
    3028       101360 :             DO idir = 1, 3
    3029        82120 :                CALL deallocate_pw(v_drho_r(idir, ispin), pw_pool)
    3030              :             END DO
    3031              :          END DO
    3032        19240 :          DEALLOCATE (v_drho, v_drho_r)
    3033              : 
    3034              :       END IF
    3035              : 
    3036        30762 :       IF (laplace_f) THEN
    3037          386 :       DO ispin = 1, nspins
    3038          386 :          CALL deallocate_pw(v_laplace(ispin), pw_pool)
    3039              :       END DO
    3040          174 :       DEALLOCATE (v_laplace)
    3041              :       END IF
    3042              : 
    3043        30762 :       IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
    3044        11576 :          CALL pw_pool%give_back_pw(tmp_g)
    3045              :       END IF
    3046              : 
    3047        30762 :       IF (ASSOCIATED(vxc_g%pw_grid) .AND. ASSOCIATED(pw_pool)) THEN
    3048        11576 :          CALL pw_pool%give_back_pw(vxc_g)
    3049              :       END IF
    3050              : 
    3051        30762 :       IF (my_compute_virial .AND. (gradient_f .OR. laplace_f)) THEN
    3052          232 :          CALL deallocate_pw(virial_pw, pw_pool)
    3053              :       END IF
    3054              : 
    3055        30762 :       CALL timestop(handle)
    3056              : 
    3057        92286 :    END SUBROUTINE xc_calc_2nd_deriv_analytical
    3058              : 
    3059              : ! **************************************************************************************************
    3060              : !> \brief allocates grids using pw_pool (if associated) or with bounds
    3061              : !> \param pw ...
    3062              : !> \param pw_pool ...
    3063              : !> \param bo ...
    3064              : ! **************************************************************************************************
    3065        89242 :    SUBROUTINE allocate_pw(pw, pw_pool, bo)
    3066              :       TYPE(pw_r3d_rs_type), INTENT(OUT)                         :: pw
    3067              :       TYPE(pw_pool_type), INTENT(IN), POINTER            :: pw_pool
    3068              :       INTEGER, DIMENSION(2, 3), INTENT(IN)               :: bo
    3069              : 
    3070        89242 :       IF (ASSOCIATED(pw_pool)) THEN
    3071        58434 :          CALL pw_pool%create_pw(pw)
    3072        58434 :          CALL pw_zero(pw)
    3073              :       ELSE
    3074       154040 :          ALLOCATE (pw%array(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
    3075     78622016 :          pw%array = 0.0_dp
    3076              :       END IF
    3077              : 
    3078        89242 :    END SUBROUTINE allocate_pw
    3079              : 
    3080              : ! **************************************************************************************************
    3081              : !> \brief deallocates grid allocated with allocate_pw
    3082              : !> \param pw ...
    3083              : !> \param pw_pool ...
    3084              : ! **************************************************************************************************
    3085        89242 :    SUBROUTINE deallocate_pw(pw, pw_pool)
    3086              :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                       :: pw
    3087              :       TYPE(pw_pool_type), INTENT(IN), POINTER            :: pw_pool
    3088              : 
    3089        89242 :       IF (ASSOCIATED(pw_pool)) THEN
    3090        58434 :          CALL pw_pool%give_back_pw(pw)
    3091              :       ELSE
    3092        30808 :          CALL pw%release()
    3093              :       END IF
    3094              : 
    3095        89242 :    END SUBROUTINE deallocate_pw
    3096              : 
    3097              : ! **************************************************************************************************
    3098              : !> \brief updates virial from first derivative w.r.t. norm_drho
    3099              : !> \param virial_pw ...
    3100              : !> \param drho ...
    3101              : !> \param drho1 ...
    3102              : !> \param deriv_data ...
    3103              : !> \param virial_xc ...
    3104              : ! **************************************************************************************************
    3105          304 :    SUBROUTINE virial_drho_drho1(virial_pw, drho, drho1, deriv_data, virial_xc)
    3106              :       TYPE(pw_r3d_rs_type), INTENT(IN)                          :: virial_pw
    3107              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN)     :: drho, drho1
    3108              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: deriv_data
    3109              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: virial_xc
    3110              : 
    3111              :       INTEGER                                            :: idir, jdir
    3112              :       REAL(KIND=dp)                                      :: tmp
    3113              : 
    3114         1216 :       DO idir = 1, 3
    3115          912 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,virial_pw,deriv_data)
    3116              :          virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*deriv_data(:, :, :)
    3117              : !$OMP END PARALLEL WORKSHARE
    3118         3952 :          DO jdir = 1, 3
    3119              :             tmp = virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
    3120         2736 :                                                               drho1(jdir)%array(:, :, :))
    3121         2736 :             virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
    3122         3648 :             virial_xc(idir, jdir) = virial_xc(idir, jdir) + tmp
    3123              :          END DO
    3124              :       END DO
    3125              : 
    3126          304 :    END SUBROUTINE virial_drho_drho1
    3127              : 
    3128              : ! **************************************************************************************************
    3129              : !> \brief Adds virial contribution from second order potential parts
    3130              : !> \param virial_pw ...
    3131              : !> \param drho ...
    3132              : !> \param v_drho ...
    3133              : !> \param virial_xc ...
    3134              : ! **************************************************************************************************
    3135          298 :    SUBROUTINE virial_drho_drho(virial_pw, drho, v_drho, virial_xc)
    3136              :       TYPE(pw_r3d_rs_type), INTENT(IN)                          :: virial_pw
    3137              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN)     :: drho
    3138              :       TYPE(pw_r3d_rs_type), INTENT(IN)                        :: v_drho
    3139              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: virial_xc
    3140              : 
    3141              :       INTEGER                                            :: idir, jdir
    3142              :       REAL(KIND=dp)                                      :: tmp
    3143              : 
    3144         1192 :       DO idir = 1, 3
    3145          894 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,idir,v_drho,virial_pw)
    3146              :          virial_pw%array(:, :, :) = drho(idir)%array(:, :, :)*v_drho%array(:, :, :)
    3147              : !$OMP END PARALLEL WORKSHARE
    3148         2980 :          DO jdir = 1, idir
    3149              :             tmp = -virial_pw%pw_grid%dvol*accurate_dot_product(virial_pw%array(:, :, :), &
    3150         1788 :                                                                drho(jdir)%array(:, :, :))
    3151         1788 :             virial_xc(jdir, idir) = virial_xc(jdir, idir) + tmp
    3152         2682 :             virial_xc(idir, jdir) = virial_xc(jdir, idir)
    3153              :          END DO
    3154              :       END DO
    3155              : 
    3156          298 :    END SUBROUTINE virial_drho_drho
    3157              : 
    3158              : ! **************************************************************************************************
    3159              : !> \brief ...
    3160              : !> \param rho_r ...
    3161              : !> \param pw_pool ...
    3162              : !> \param virial_xc ...
    3163              : !> \param deriv_data ...
    3164              : ! **************************************************************************************************
    3165          150 :    SUBROUTINE virial_laplace(rho_r, pw_pool, virial_xc, deriv_data)
    3166              :       TYPE(pw_r3d_rs_type), TARGET :: rho_r
    3167              :       TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
    3168              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: virial_xc
    3169              :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: deriv_data
    3170              : 
    3171              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'virial_laplace'
    3172              : 
    3173              :       INTEGER :: handle, idir, jdir
    3174              :       TYPE(pw_r3d_rs_type), POINTER :: virial_pw
    3175              :       TYPE(pw_c1d_gs_type), POINTER :: tmp_g, rho_g
    3176              :       INTEGER, DIMENSION(3) :: my_deriv
    3177              : 
    3178          150 :       CALL timeset(routineN, handle)
    3179              : 
    3180          150 :       NULLIFY (virial_pw, tmp_g, rho_g)
    3181          150 :       ALLOCATE (virial_pw, tmp_g, rho_g)
    3182          150 :       CALL pw_pool%create_pw(virial_pw)
    3183          150 :       CALL pw_pool%create_pw(tmp_g)
    3184          150 :       CALL pw_pool%create_pw(rho_g)
    3185          150 :       CALL pw_zero(virial_pw)
    3186          150 :       CALL pw_transfer(rho_r, rho_g)
    3187          600 :       DO idir = 1, 3
    3188         1500 :          DO jdir = idir, 3
    3189          900 :             CALL pw_copy(rho_g, tmp_g)
    3190              : 
    3191          900 :             my_deriv = 0
    3192          900 :             my_deriv(idir) = 1
    3193          900 :             my_deriv(jdir) = my_deriv(jdir) + 1
    3194              : 
    3195          900 :             CALL pw_derive(tmp_g, my_deriv)
    3196          900 :             CALL pw_transfer(tmp_g, virial_pw)
    3197              :             virial_xc(idir, jdir) = virial_xc(idir, jdir) - 2.0_dp*virial_pw%pw_grid%dvol* &
    3198              :                                     accurate_dot_product(virial_pw%array(:, :, :), &
    3199          900 :                                                          deriv_data(:, :, :))
    3200         1350 :             virial_xc(jdir, idir) = virial_xc(idir, jdir)
    3201              :          END DO
    3202              :       END DO
    3203          150 :       CALL pw_pool%give_back_pw(virial_pw)
    3204          150 :       CALL pw_pool%give_back_pw(tmp_g)
    3205          150 :       CALL pw_pool%give_back_pw(rho_g)
    3206          150 :       DEALLOCATE (virial_pw, tmp_g, rho_g)
    3207              : 
    3208          150 :       CALL timestop(handle)
    3209              : 
    3210          150 :    END SUBROUTINE virial_laplace
    3211              : 
    3212              : ! **************************************************************************************************
    3213              : !> \brief Prepare objects for the calculation of the 2nd derivatives of the density functional.
    3214              : !>      The calculation must then be performed with xc_calc_2nd_deriv.
    3215              : !> \param deriv_set object containing the XC derivatives (out)
    3216              : !> \param rho_set object that will contain the density at which the
    3217              : !>        derivatives were calculated
    3218              : !> \param rho_r the place where you evaluate the derivative
    3219              : !> \param pw_pool the pool for the grids
    3220              : !> \param xc_section which functional should be used and how to calculate it
    3221              : !> \param tau_r kinetic energy density in real space
    3222              : ! **************************************************************************************************
    3223        12134 :    SUBROUTINE xc_prep_2nd_deriv(deriv_set, &
    3224              :                                 rho_set, rho_r, pw_pool, xc_section, tau_r)
    3225              : 
    3226              :       TYPE(xc_derivative_set_type)                       :: deriv_set
    3227              :       TYPE(xc_rho_set_type)                              :: rho_set
    3228              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER               :: rho_r
    3229              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
    3230              :       TYPE(section_vals_type), POINTER                   :: xc_section
    3231              :       TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, POINTER     :: tau_r
    3232              : 
    3233              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xc_prep_2nd_deriv'
    3234              : 
    3235              :       INTEGER                                            :: handle, nspins
    3236              :       LOGICAL                                            :: lsd
    3237        12134 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER               :: rho_g
    3238        12134 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau
    3239              : 
    3240        12134 :       CALL timeset(routineN, handle)
    3241              : 
    3242        12134 :       CPASSERT(ASSOCIATED(xc_section))
    3243        12134 :       CPASSERT(ASSOCIATED(pw_pool))
    3244              : 
    3245        12134 :       nspins = SIZE(rho_r)
    3246        12134 :       lsd = (nspins /= 1)
    3247              : 
    3248        12134 :       NULLIFY (rho_g, tau)
    3249        12134 :       IF (PRESENT(tau_r)) &
    3250        11472 :          tau => tau_r
    3251              : 
    3252        12134 :       IF (section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")) THEN
    3253              :          CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 2, &
    3254              :                                          rho_r, rho_g, tau, xc_section, pw_pool, &
    3255        12050 :                                          calc_potential=.TRUE.)
    3256              :       ELSE
    3257              :          CALL xc_rho_set_and_dset_create(rho_set, deriv_set, 1, &
    3258              :                                          rho_r, rho_g, tau, xc_section, pw_pool, &
    3259           84 :                                          calc_potential=.TRUE.)
    3260              :       END IF
    3261              : 
    3262        12134 :       CALL timestop(handle)
    3263              : 
    3264        12134 :    END SUBROUTINE xc_prep_2nd_deriv
    3265              : 
    3266              : ! **************************************************************************************************
    3267              : !> \brief divides derivatives from deriv_set by norm_drho
    3268              : !> \param deriv_set ...
    3269              : !> \param rho_set ...
    3270              : !> \param lsd ...
    3271              : ! **************************************************************************************************
    3272       144081 :    SUBROUTINE divide_by_norm_drho(deriv_set, rho_set, lsd)
    3273              : 
    3274              :       TYPE(xc_derivative_set_type), INTENT(INOUT)        :: deriv_set
    3275              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
    3276              :       LOGICAL, INTENT(IN)                                :: lsd
    3277              : 
    3278       144081 :       INTEGER, DIMENSION(:), POINTER                     :: split_desc
    3279              :       INTEGER                                            :: idesc
    3280              :       INTEGER, DIMENSION(2, 3)                           :: bo
    3281              :       REAL(KIND=dp)                                      :: drho_cutoff
    3282       144081 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: norm_drho, norm_drhoa, norm_drhob
    3283      1728972 :       TYPE(cp_3d_r_cp_type), DIMENSION(3)                 :: drho, drhoa, drhob
    3284              :       TYPE(cp_sll_xc_deriv_type), POINTER                :: pos
    3285              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att
    3286              : 
    3287              : ! check for unknown derivatives and divide by norm_drho where necessary
    3288              : 
    3289      1440810 :       bo = rho_set%local_bounds
    3290              :       CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, norm_drho=norm_drho, &
    3291              :                           norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
    3292       144081 :                           drho=drho, drhoa=drhoa, drhob=drhob, can_return_null=.TRUE.)
    3293              : 
    3294       144081 :       pos => deriv_set%derivs
    3295       611467 :       DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
    3296       467386 :          CALL xc_derivative_get(deriv_att, split_desc=split_desc)
    3297       995792 :          DO idesc = 1, SIZE(split_desc)
    3298       467386 :             SELECT CASE (split_desc(idesc))
    3299              :             CASE (deriv_norm_drho)
    3300       116801 :                IF (ASSOCIATED(norm_drho)) THEN
    3301       116801 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drho,drho_cutoff)
    3302              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3303              :                                                   MAX(norm_drho(:, :, :), drho_cutoff)
    3304              : !$OMP END PARALLEL WORKSHARE
    3305            0 :                ELSE IF (ASSOCIATED(drho(1)%array)) THEN
    3306            0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drho,drho_cutoff)
    3307              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3308              :                                                   MAX(SQRT(drho(1)%array(:, :, :)**2 + &
    3309              :                                                            drho(2)%array(:, :, :)**2 + &
    3310              :                                                            drho(3)%array(:, :, :)**2), drho_cutoff)
    3311              : !$OMP END PARALLEL WORKSHARE
    3312            0 :                ELSE IF (ASSOCIATED(drhoa(1)%array) .AND. ASSOCIATED(drhob(1)%array)) THEN
    3313            0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drhob,drho_cutoff)
    3314              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3315              :                                                   MAX(SQRT((drhoa(1)%array(:, :, :) + drhob(1)%array(:, :, :))**2 + &
    3316              :                                                            (drhoa(2)%array(:, :, :) + drhob(2)%array(:, :, :))**2 + &
    3317              :                                                            (drhoa(3)%array(:, :, :) + drhob(3)%array(:, :, :))**2), drho_cutoff)
    3318              : !$OMP END PARALLEL WORKSHARE
    3319              :                ELSE
    3320            0 :                   CPABORT("Normalization of derivative requires any of norm_drho, drho or drhoa+drhob!")
    3321              :                END IF
    3322              :             CASE (deriv_norm_drhoa)
    3323        19264 :                IF (ASSOCIATED(norm_drhoa)) THEN
    3324        19264 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhoa,drho_cutoff)
    3325              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3326              :                                                   MAX(norm_drhoa(:, :, :), drho_cutoff)
    3327              : !$OMP END PARALLEL WORKSHARE
    3328            0 :                ELSE IF (ASSOCIATED(drhoa(1)%array)) THEN
    3329            0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhoa,drho_cutoff)
    3330              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3331              :                                                   MAX(SQRT(drhoa(1)%array(:, :, :)**2 + &
    3332              :                                                            drhoa(2)%array(:, :, :)**2 + &
    3333              :                                                            drhoa(3)%array(:, :, :)**2), drho_cutoff)
    3334              : !$OMP END PARALLEL WORKSHARE
    3335              :                ELSE
    3336            0 :                   CPABORT("Normalization of derivative requires any of norm_drhoa or drhoa!")
    3337              :                END IF
    3338              :             CASE (deriv_norm_drhob)
    3339        19260 :                IF (ASSOCIATED(norm_drhob)) THEN
    3340        19260 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,norm_drhob,drho_cutoff)
    3341              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3342              :                                                   MAX(norm_drhob(:, :, :), drho_cutoff)
    3343              : !$OMP END PARALLEL WORKSHARE
    3344            0 :                ELSE IF (ASSOCIATED(drhob(1)%array)) THEN
    3345            0 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(deriv_att,drhob,drho_cutoff)
    3346              :                   deriv_att%deriv_data(:, :, :) = deriv_att%deriv_data(:, :, :)/ &
    3347              :                                                   MAX(SQRT(drhob(1)%array(:, :, :)**2 + &
    3348              :                                                            drhob(2)%array(:, :, :)**2 + &
    3349              :                                                            drhob(3)%array(:, :, :)**2), drho_cutoff)
    3350              : !$OMP END PARALLEL WORKSHARE
    3351              :                ELSE
    3352            0 :                   CPABORT("Normalization of derivative requires any of norm_drhob or drhob!")
    3353              :                END IF
    3354              :             CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
    3355       162288 :                IF (lsd) &
    3356            0 :                   CPABORT(TRIM(id_to_desc(split_desc(idesc)))//" not handled in lsd!'")
    3357              :             CASE (deriv_rhoa, deriv_rhob, deriv_tau_a, deriv_tau_b, deriv_laplace_rhoa, deriv_laplace_rhob)
    3358              :             CASE default
    3359       384325 :                CPABORT("Unknown derivative id")
    3360              :             END SELECT
    3361              :          END DO
    3362              :       END DO
    3363              : 
    3364       144081 :    END SUBROUTINE divide_by_norm_drho
    3365              : 
    3366              : ! **************************************************************************************************
    3367              : !> \brief allocates and calculates drho from given spin densities drhoa, drhob
    3368              : !> \param drho ...
    3369              : !> \param drhoa ...
    3370              : !> \param drhob ...
    3371              : ! **************************************************************************************************
    3372        16080 :    SUBROUTINE calc_drho_from_ab(drho, drhoa, drhob)
    3373              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(OUT)    :: drho
    3374              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN)     :: drhoa, drhob
    3375              : 
    3376              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calc_drho_from_ab'
    3377              : 
    3378              :       INTEGER                                            :: handle, idir
    3379              : 
    3380         4020 :       CALL timeset(routineN, handle)
    3381              : 
    3382        16080 :       DO idir = 1, 3
    3383        12060 :          NULLIFY (drho(idir)%array)
    3384              :          ALLOCATE (drho(idir)%array(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
    3385              :                                     LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
    3386       144720 :                                     LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
    3387        16080 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(drho,drhoa,drhob,idir)
    3388              :          drho(idir)%array(:, :, :) = drhoa(idir)%array(:, :, :) + drhob(idir)%array(:, :, :)
    3389              : !$OMP END PARALLEL WORKSHARE
    3390              :       END DO
    3391              : 
    3392         4020 :       CALL timestop(handle)
    3393              : 
    3394         4020 :    END SUBROUTINE calc_drho_from_ab
    3395              : 
    3396              : ! **************************************************************************************************
    3397              : !> \brief allocates and calculates dot products of two density gradients
    3398              : !> \param dr1dr ...
    3399              : !> \param drho ...
    3400              : !> \param drho1 ...
    3401              : ! **************************************************************************************************
    3402        22586 :    SUBROUTINE prepare_dr1dr(dr1dr, drho, drho1)
    3403              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
    3404              :          INTENT(OUT)                                     :: dr1dr
    3405              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN)     :: drho, drho1
    3406              : 
    3407              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'prepare_dr1dr'
    3408              : 
    3409              :       INTEGER                                            :: handle, idir
    3410              : 
    3411        22586 :       CALL timeset(routineN, handle)
    3412              : 
    3413            0 :       ALLOCATE (dr1dr(LBOUND(drho(1)%array, 1):UBOUND(drho(1)%array, 1), &
    3414              :                       LBOUND(drho(1)%array, 2):UBOUND(drho(1)%array, 2), &
    3415       271032 :                       LBOUND(drho(1)%array, 3):UBOUND(drho(1)%array, 3)))
    3416              : 
    3417        22586 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1)
    3418              :       dr1dr(:, :, :) = drho(1)%array(:, :, :)*drho1(1)%array(:, :, :)
    3419              : !$OMP END PARALLEL WORKSHARE
    3420        67758 :       DO idir = 2, 3
    3421        67758 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(dr1dr,drho,drho1,idir)
    3422              :          dr1dr(:, :, :) = dr1dr(:, :, :) + drho(idir)%array(:, :, :)*drho1(idir)%array(:, :, :)
    3423              : !$OMP END PARALLEL WORKSHARE
    3424              :       END DO
    3425              : 
    3426        22586 :       CALL timestop(handle)
    3427              : 
    3428        22586 :    END SUBROUTINE prepare_dr1dr
    3429              : 
    3430              : ! **************************************************************************************************
    3431              : !> \brief allocates and calculates dot product of two densities for triplets
    3432              : !> \param dr1dr ...
    3433              : !> \param drhoa ...
    3434              : !> \param drhob ...
    3435              : !> \param drho1a ...
    3436              : !> \param drho1b ...
    3437              : !> \param fac ...
    3438              : ! **************************************************************************************************
    3439          710 :    SUBROUTINE prepare_dr1dr_ab(dr1dr, drhoa, drhob, drho1a, drho1b, fac)
    3440              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
    3441              :          INTENT(OUT)                                     :: dr1dr
    3442              :       TYPE(cp_3d_r_cp_type), DIMENSION(3), INTENT(IN)    :: drhoa, drhob, drho1a, drho1b
    3443              :       REAL(KIND=dp), INTENT(IN)                          :: fac
    3444              : 
    3445              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'prepare_dr1dr_ab'
    3446              : 
    3447              :       INTEGER                                            :: handle, idir
    3448              : 
    3449          710 :       CALL timeset(routineN, handle)
    3450              : 
    3451            0 :       ALLOCATE (dr1dr(LBOUND(drhoa(1)%array, 1):UBOUND(drhoa(1)%array, 1), &
    3452              :                       LBOUND(drhoa(1)%array, 2):UBOUND(drhoa(1)%array, 2), &
    3453         8520 :                       LBOUND(drhoa(1)%array, 3):UBOUND(drhoa(1)%array, 3)))
    3454              : 
    3455          710 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob)
    3456              :       dr1dr(:, :, :) = drhoa(1)%array(:, :, :)*(drho1a(1)%array(:, :, :) + &
    3457              :                                                 fac*drho1b(1)%array(:, :, :)) + &
    3458              :                        drhob(1)%array(:, :, :)*(fac*drho1a(1)%array(:, :, :) + &
    3459              :                                                 drho1b(1)%array(:, :, :))
    3460              : !$OMP END PARALLEL WORKSHARE
    3461         2130 :       DO idir = 2, 3
    3462         2130 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(fac,dr1dr,drho1a,drho1b,drhoa,drhob,idir)
    3463              :          dr1dr(:, :, :) = dr1dr(:, :, :) + &
    3464              :                           drhoa(idir)%array(:, :, :)*(drho1a(idir)%array(:, :, :) + &
    3465              :                                                       fac*drho1b(idir)%array(:, :, :)) + &
    3466              :                           drhob(idir)%array(:, :, :)*(fac*drho1a(idir)%array(:, :, :) + &
    3467              :                                                       drho1b(idir)%array(:, :, :))
    3468              : !$OMP END PARALLEL WORKSHARE
    3469              :       END DO
    3470              : 
    3471          710 :       CALL timestop(handle)
    3472              : 
    3473          710 :    END SUBROUTINE prepare_dr1dr_ab
    3474              : 
    3475              : ! **************************************************************************************************
    3476              : !> \brief checks for gradients
    3477              : !> \param deriv_set ...
    3478              : !> \param lsd ...
    3479              : !> \param gradient_f ...
    3480              : !> \param tau_f ...
    3481              : !> \param laplace_f ...
    3482              : ! **************************************************************************************************
    3483       142037 :    SUBROUTINE check_for_derivatives(deriv_set, lsd, rho_f, gradient_f, tau_f, laplace_f)
    3484              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
    3485              :       LOGICAL, INTENT(IN)                                :: lsd
    3486              :       LOGICAL, INTENT(OUT)                               :: rho_f, gradient_f, tau_f, laplace_f
    3487              : 
    3488              :       CHARACTER(len=*), PARAMETER :: routineN = 'check_for_derivatives'
    3489              : 
    3490              :       INTEGER                                            :: handle, iorder, order
    3491       142037 :       INTEGER, DIMENSION(:), POINTER                     :: split_desc
    3492              :       TYPE(cp_sll_xc_deriv_type), POINTER                :: pos
    3493              :       TYPE(xc_derivative_type), POINTER                  :: deriv_att
    3494              : 
    3495       142037 :       CALL timeset(routineN, handle)
    3496              : 
    3497       142037 :       rho_f = .FALSE.
    3498       142037 :       gradient_f = .FALSE.
    3499       142037 :       tau_f = .FALSE.
    3500       142037 :       laplace_f = .FALSE.
    3501              :       ! check for unknown derivatives
    3502       142037 :       pos => deriv_set%derivs
    3503       645375 :       DO WHILE (cp_sll_xc_deriv_next(pos, el_att=deriv_att))
    3504              :          CALL xc_derivative_get(deriv_att, order=order, &
    3505       503338 :                                 split_desc=split_desc)
    3506       645375 :          IF (lsd) THEN
    3507       301081 :             DO iorder = 1, size(split_desc)
    3508       150257 :                SELECT CASE (split_desc(iorder))
    3509              :                CASE (deriv_rhoa, deriv_rhob)
    3510        77920 :                   rho_f = .TRUE.
    3511              :                CASE (deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob)
    3512        68832 :                   gradient_f = .TRUE.
    3513              :                CASE (deriv_tau_a, deriv_tau_b)
    3514         2760 :                   tau_f = .TRUE.
    3515              :                CASE (deriv_laplace_rhoa, deriv_laplace_rhob)
    3516         1312 :                   laplace_f = .TRUE.
    3517              :                CASE (deriv_rho, deriv_tau, deriv_laplace_rho)
    3518            0 :                   CPABORT("Derivative not handled in lsd!")
    3519              :                CASE default
    3520       150824 :                   CPABORT("Unknown derivative id")
    3521              :                END SELECT
    3522              :             END DO
    3523              :          ELSE
    3524       653818 :             DO iorder = 1, size(split_desc)
    3525       353081 :                SELECT CASE (split_desc(iorder))
    3526              :                CASE (deriv_rho)
    3527       177956 :                   rho_f = .TRUE.
    3528              :                CASE (deriv_tau)
    3529         4948 :                   tau_f = .TRUE.
    3530              :                CASE (deriv_norm_drho)
    3531       116395 :                   gradient_f = .TRUE.
    3532              :                CASE (deriv_laplace_rho)
    3533         1438 :                   laplace_f = .TRUE.
    3534              :                CASE default
    3535       300737 :                   CPABORT("Unknown derivative id")
    3536              :                END SELECT
    3537              :             END DO
    3538              :          END IF
    3539              :       END DO
    3540              : 
    3541       142037 :       CALL timestop(handle)
    3542              : 
    3543       142037 :    END SUBROUTINE check_for_derivatives
    3544              : 
    3545              : END MODULE xc
        

Generated by: LCOV version 2.0-1