LCOV - code coverage report
Current view: top level - src/xc - xc_rho_set_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b279b6b) Lines: 353 358 98.6 %
Date: 2024-04-24 07:13:09 Functions: 7 8 87.5 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief contains the structure
      10             : !> \par History
      11             : !>      11.2003 created [fawzi]
      12             : !> \author fawzi
      13             : ! **************************************************************************************************
      14             : MODULE xc_rho_set_types
      15             :    USE cp_array_utils, ONLY: cp_3d_r_cp_type
      16             :    USE kinds, ONLY: dp
      17             :    USE pw_grid_types, ONLY: pw_grid_type
      18             :    USE pw_methods, ONLY: pw_copy, &
      19             :                          pw_transfer
      20             :    USE pw_pool_types, ONLY: &
      21             :       pw_pool_type
      22             :    USE pw_spline_utils, ONLY: pw_spline_scale_deriv
      23             :    USE pw_types, ONLY: &
      24             :       pw_c1d_gs_type, &
      25             :       pw_r3d_rs_type
      26             :    USE xc_input_constants, ONLY: xc_deriv_pw, &
      27             :                                  xc_deriv_spline2, &
      28             :                                  xc_deriv_spline3, &
      29             :                                  xc_rho_no_smooth
      30             :    USE xc_rho_cflags_types, ONLY: xc_rho_cflags_equal, &
      31             :                                   xc_rho_cflags_setall, &
      32             :                                   xc_rho_cflags_type
      33             :    USE xc_util, ONLY: xc_pw_gradient, &
      34             :                       xc_pw_laplace, &
      35             :                       xc_pw_smooth
      36             : #include "../base/base_uses.f90"
      37             : 
      38             :    IMPLICIT NONE
      39             :    PRIVATE
      40             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      41             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_rho_set_types'
      42             : 
      43             :    PUBLIC :: xc_rho_set_type
      44             :    PUBLIC :: xc_rho_set_create, xc_rho_set_release, &
      45             :              xc_rho_set_update, xc_rho_set_get, xc_rho_set_recover_pw
      46             : 
      47             : ! **************************************************************************************************
      48             : !> \brief represent a density, with all the representation and data needed
      49             : !>      to perform a functional evaluation
      50             : !> \param local_bounds the part of the 3d array on which the functional is
      51             : !>        computed
      52             : !> \param owns which components are owned by this structure (and should
      53             : !>        be deallocated
      54             : !> \param has which components are present and up to date
      55             : !> \param rho the density
      56             : !> \param drho the gradient of the density (x,y and z direction)
      57             : !> \param norm_drho the norm of the gradient of the density
      58             : !> \param rhoa , rhob: spin alpha and beta parts of the density in the LSD case
      59             : !> \param drhoa , drhob: gradient of the spin alpha and beta parts of the density
      60             : !>        in the LSD case (x,y and z direction)
      61             : !> \param norm_drhoa , norm_drhob: norm of the gradient of rhoa and rhob
      62             : !> \param rho_ 1_3: rho^(1.0_dp/3.0_dp)
      63             : !> \param rhoa_ 1_3, rhob_1_3: rhoa^(1.0_dp/3.0_dp), rhob^(1.0_dp/3.0_dp)
      64             : !> \param tau the kinetic (KohnSham) part of rho
      65             : !> \param tau_a the kinetic (KohnSham) part of rhoa
      66             : !> \param tau_b the kinetic (KohnSham) part of rhob
      67             : !> \note
      68             : !>      the use of 3d arrays is the result of trying to use only basic
      69             : !>      types (to be generic and independent from the method), and
      70             : !>      avoiding copies using the actual structure.
      71             : !>      only the part defined by local bounds is guaranteed to be present,
      72             : !>      and it is the only meaningful part.
      73             : !> \par History
      74             : !>      11.2003 created [fawzi & thomas]
      75             : !>      12.2008 added laplace parts [mguidon]
      76             : !> \author fawzi & thomas
      77             : ! **************************************************************************************************
      78             :    TYPE xc_rho_set_type
      79             :       INTEGER, DIMENSION(2, 3) :: local_bounds = -1
      80             :       REAL(kind=dp) :: rho_cutoff = EPSILON(0.0_dp), drho_cutoff = EPSILON(0.0_dp), tau_cutoff = EPSILON(0.0_dp)
      81             :       TYPE(xc_rho_cflags_type) :: owns = xc_rho_cflags_type(), has = xc_rho_cflags_type()
      82             :       ! for spin restricted systems
      83             :       REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rho => NULL()
      84             :       TYPE(cp_3d_r_cp_type), DIMENSION(3)         :: drho = cp_3d_r_cp_type()
      85             :       REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: norm_drho => NULL()
      86             :       REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rho_1_3 => NULL()
      87             :       REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: tau => NULL()
      88             :       ! for UNrestricted systems
      89             :       REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rhoa => NULL(), rhob => NULL()
      90             :       TYPE(cp_3d_r_cp_type), DIMENSION(3)         :: drhoa = cp_3d_r_cp_type(), &
      91             :                                                      drhob = cp_3d_r_cp_type()
      92             :       REAL(KIND=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: norm_drhoa => NULL(), norm_drhob => NULL()
      93             :       REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: rhoa_1_3 => NULL(), rhob_1_3 => NULL()
      94             :       REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: tau_a => NULL(), tau_b => NULL()
      95             :       REAL(kind=dp), DIMENSION(:, :, :), CONTIGUOUS, POINTER :: laplace_rho => NULL(), laplace_rhoa => NULL(), &
      96             :                                                                 laplace_rhob => NULL()
      97             :    END TYPE xc_rho_set_type
      98             : 
      99             : CONTAINS
     100             : 
     101             : ! **************************************************************************************************
     102             : !> \brief allocates and does (minimal) initialization of a rho_set
     103             : !> \param rho_set the structure to allocate
     104             : !> \param local_bounds ...
     105             : !> \param rho_cutoff ...
     106             : !> \param drho_cutoff ...
     107             : !> \param tau_cutoff ...
     108             : ! **************************************************************************************************
     109      248727 :    SUBROUTINE xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, &
     110             :                                 tau_cutoff)
     111             :       TYPE(xc_rho_set_type), INTENT(INOUT)               :: rho_set
     112             :       INTEGER, DIMENSION(2, 3), INTENT(in)               :: local_bounds
     113             :       REAL(kind=dp), INTENT(in), OPTIONAL                :: rho_cutoff, drho_cutoff, tau_cutoff
     114             : 
     115      248727 :       IF (PRESENT(rho_cutoff)) rho_set%rho_cutoff = rho_cutoff
     116      248727 :       IF (PRESENT(drho_cutoff)) rho_set%drho_cutoff = drho_cutoff
     117      248727 :       IF (PRESENT(tau_cutoff)) rho_set%tau_cutoff = tau_cutoff
     118     2487270 :       rho_set%local_bounds = local_bounds
     119      248727 :       CALL xc_rho_cflags_setall(rho_set%owns, .TRUE.)
     120      248727 :       CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
     121      248727 :    END SUBROUTINE xc_rho_set_create
     122             : 
     123             : ! **************************************************************************************************
     124             : !> \brief releases the given rho_set
     125             : !> \param rho_set the structure to release
     126             : !> \param pw_pool the plae where to give back the arrays
     127             : ! **************************************************************************************************
     128      248727 :    SUBROUTINE xc_rho_set_release(rho_set, pw_pool)
     129             :       TYPE(xc_rho_set_type), INTENT(INOUT)               :: rho_set
     130             :       TYPE(pw_pool_type), OPTIONAL, POINTER              :: pw_pool
     131             : 
     132             :       INTEGER                                            :: i
     133             : 
     134      248727 :       IF (PRESENT(pw_pool)) THEN
     135      119067 :          IF (ASSOCIATED(pw_pool)) THEN
     136      119067 :             CALL xc_rho_set_clean(rho_set, pw_pool)
     137             :          ELSE
     138           0 :             CPABORT("pw_pool must be associated")
     139             :          END IF
     140             :       END IF
     141             : 
     142      994908 :       rho_set%local_bounds(1, :) = -HUGE(0) ! we want to crash...
     143      994908 :       rho_set%local_bounds(1, :) = HUGE(0)
     144      248727 :       IF (rho_set%owns%rho .AND. ASSOCIATED(rho_set%rho)) THEN
     145      112305 :          DEALLOCATE (rho_set%rho)
     146             :       ELSE
     147      136422 :          NULLIFY (rho_set%rho)
     148             :       END IF
     149      248727 :       IF (rho_set%owns%rho_spin) THEN
     150      108758 :          IF (ASSOCIATED(rho_set%rhoa)) THEN
     151       17277 :             DEALLOCATE (rho_set%rhoa)
     152             :          END IF
     153      108758 :          IF (ASSOCIATED(rho_set%rhob)) THEN
     154       17277 :             DEALLOCATE (rho_set%rhob)
     155             :          END IF
     156             :       ELSE
     157      139969 :          NULLIFY (rho_set%rhoa, rho_set%rhob)
     158             :       END IF
     159      248727 :       IF (rho_set%owns%rho_1_3 .AND. ASSOCIATED(rho_set%rho_1_3)) THEN
     160        5345 :          DEALLOCATE (rho_set%rho_1_3)
     161             :       ELSE
     162      243382 :          NULLIFY (rho_set%rho_1_3)
     163             :       END IF
     164      248727 :       IF (rho_set%owns%rho_spin) THEN
     165      108758 :          IF (ASSOCIATED(rho_set%rhoa_1_3)) THEN
     166        2456 :             DEALLOCATE (rho_set%rhoa_1_3)
     167             :          END IF
     168      108758 :          IF (ASSOCIATED(rho_set%rhob_1_3)) THEN
     169        2456 :             DEALLOCATE (rho_set%rhob_1_3)
     170             :          END IF
     171             :       ELSE
     172      139969 :          NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
     173             :       END IF
     174      248727 :       IF (rho_set%owns%drho) THEN
     175      471168 :          DO i = 1, 3
     176      471168 :             IF (ASSOCIATED(rho_set%drho(i)%array)) THEN
     177      173484 :                DEALLOCATE (rho_set%drho(i)%array)
     178             :             END IF
     179             :          END DO
     180             :       ELSE
     181      523740 :          DO i = 1, 3
     182      523740 :             NULLIFY (rho_set%drho(i)%array)
     183             :          END DO
     184             :       END IF
     185      248727 :       IF (rho_set%owns%drho_spin) THEN
     186      428528 :          DO i = 1, 3
     187      321396 :             IF (ASSOCIATED(rho_set%drhoa(i)%array)) THEN
     188       28254 :                DEALLOCATE (rho_set%drhoa(i)%array)
     189             :             END IF
     190      428528 :             IF (ASSOCIATED(rho_set%drhob(i)%array)) THEN
     191       28254 :                DEALLOCATE (rho_set%drhob(i)%array)
     192             :             END IF
     193             :          END DO
     194             :       ELSE
     195      566380 :          DO i = 1, 3
     196      566380 :             NULLIFY (rho_set%drhoa(i)%array, rho_set%drhob(i)%array)
     197             :          END DO
     198             :       END IF
     199      248727 :       IF (rho_set%owns%laplace_rho .AND. ASSOCIATED(rho_set%laplace_rho)) THEN
     200         190 :          DEALLOCATE (rho_set%laplace_rho)
     201             :       ELSE
     202      248537 :          NULLIFY (rho_set%laplace_rho)
     203             :       END IF
     204             : 
     205      248727 :       IF (rho_set%owns%norm_drho .AND. ASSOCIATED(rho_set%norm_drho)) THEN
     206       86786 :          DEALLOCATE (rho_set%norm_drho)
     207             :       ELSE
     208      161941 :          NULLIFY (rho_set%norm_drho)
     209             :       END IF
     210      248727 :       IF (rho_set%owns%laplace_rho_spin) THEN
     211      105652 :          IF (ASSOCIATED(rho_set%laplace_rhoa)) THEN
     212          36 :             DEALLOCATE (rho_set%laplace_rhoa)
     213             :          END IF
     214      105652 :          IF (ASSOCIATED(rho_set%laplace_rhob)) THEN
     215          36 :             DEALLOCATE (rho_set%laplace_rhob)
     216             :          END IF
     217             :       ELSE
     218      143075 :          NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
     219             :       END IF
     220             : 
     221      248727 :       IF (rho_set%owns%norm_drho_spin) THEN
     222      107132 :          IF (ASSOCIATED(rho_set%norm_drhoa)) THEN
     223       10297 :             DEALLOCATE (rho_set%norm_drhoa)
     224             :          END IF
     225      107132 :          IF (ASSOCIATED(rho_set%norm_drhob)) THEN
     226       10297 :             DEALLOCATE (rho_set%norm_drhob)
     227             :          END IF
     228             :       ELSE
     229      141595 :          NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
     230             :       END IF
     231      248727 :       IF (rho_set%owns%tau .AND. ASSOCIATED(rho_set%tau)) THEN
     232         996 :          DEALLOCATE (rho_set%tau)
     233             :       ELSE
     234      247731 :          NULLIFY (rho_set%tau)
     235             :       END IF
     236      248727 :       IF (rho_set%owns%tau_spin) THEN
     237      105616 :          IF (ASSOCIATED(rho_set%tau_a)) THEN
     238           2 :             DEALLOCATE (rho_set%tau_a)
     239             :          END IF
     240      105616 :          IF (ASSOCIATED(rho_set%tau_b)) THEN
     241           2 :             DEALLOCATE (rho_set%tau_b)
     242             :          END IF
     243             :       ELSE
     244      143111 :          NULLIFY (rho_set%tau_a, rho_set%tau_b)
     245             :       END IF
     246      248727 :    END SUBROUTINE xc_rho_set_release
     247             : 
     248             : ! **************************************************************************************************
     249             : !> \brief returns the various attributes of rho_set
     250             : !> \param rho_set the object you want info about
     251             : !> \param can_return_null if true the object returned can be null,
     252             : !>        if false (the default) it stops with an error if a requested
     253             : !>        component is not associated
     254             : !> \param rho ...
     255             : !> \param drho ...
     256             : !> \param norm_drho ...
     257             : !> \param rhoa ...
     258             : !> \param rhob ...
     259             : !> \param norm_drhoa ...
     260             : !> \param norm_drhob ...
     261             : !> \param rho_1_3 ...
     262             : !> \param rhoa_1_3 ...
     263             : !> \param rhob_1_3 ...
     264             : !> \param laplace_rho ...
     265             : !> \param laplace_rhoa ...
     266             : !> \param laplace_rhob ...
     267             : !> \param drhoa ...
     268             : !> \param drhob ...
     269             : !> \param rho_cutoff ...
     270             : !> \param drho_cutoff ...
     271             : !> \param tau_cutoff ...
     272             : !> \param tau ...
     273             : !> \param tau_a ...
     274             : !> \param tau_b ...
     275             : !> \param local_bounds ...
     276             : ! **************************************************************************************************
     277     1004586 :    SUBROUTINE xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, &
     278             :                              rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, &
     279             :                              rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, &
     280             :                              drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
     281             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     282             :       LOGICAL, INTENT(in), OPTIONAL                      :: can_return_null
     283             :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     284             :          POINTER                                         :: rho
     285             :       TYPE(cp_3d_r_cp_type), DIMENSION(3), OPTIONAL       :: drho
     286             :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     287             :          POINTER                                         :: norm_drho, rhoa, rhob, norm_drhoa, &
     288             :                                                             norm_drhob, rho_1_3, rhoa_1_3, &
     289             :                                                             rhob_1_3, laplace_rho, laplace_rhoa, &
     290             :                                                             laplace_rhob
     291             :       TYPE(cp_3d_r_cp_type), DIMENSION(3), OPTIONAL       :: drhoa, drhob
     292             :       REAL(kind=dp), INTENT(out), OPTIONAL               :: rho_cutoff, drho_cutoff, tau_cutoff
     293             :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     294             :          POINTER                                         :: tau, tau_a, tau_b
     295             :       INTEGER, DIMENSION(2, 3), INTENT(OUT), OPTIONAL    :: local_bounds
     296             : 
     297             :       INTEGER                                            :: i
     298             :       LOGICAL                                            :: my_can_return_null
     299             : 
     300     1004586 :       my_can_return_null = .FALSE.
     301     1004586 :       IF (PRESENT(can_return_null)) my_can_return_null = can_return_null
     302             : 
     303     1004586 :       IF (PRESENT(rho)) THEN
     304      260853 :          rho => rho_set%rho
     305      260853 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(rho))
     306             :       END IF
     307     1004586 :       IF (PRESENT(drho)) THEN
     308      695164 :          DO i = 1, 3
     309      521373 :             drho(i)%array => rho_set%drho(i)%array
     310      695164 :             CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drho(i)%array))
     311             :          END DO
     312             :       END IF
     313     1004586 :       IF (PRESENT(norm_drho)) THEN
     314      361602 :          norm_drho => rho_set%norm_drho
     315      361602 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drho))
     316             :       END IF
     317     1004586 :       IF (PRESENT(laplace_rho)) THEN
     318       12228 :          laplace_rho => rho_set%laplace_rho
     319       12228 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rho))
     320             :       END IF
     321     1004586 :       IF (PRESENT(rhoa)) THEN
     322      143681 :          rhoa => rho_set%rhoa
     323      143681 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa))
     324             :       END IF
     325     1004586 :       IF (PRESENT(rhob)) THEN
     326      143681 :          rhob => rho_set%rhob
     327      143681 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob))
     328             :       END IF
     329     1004586 :       IF (PRESENT(drhoa)) THEN
     330      572028 :          DO i = 1, 3
     331      429021 :             drhoa(i)%array => rho_set%drhoa(i)%array
     332      572028 :             CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhoa(i)%array))
     333             :          END DO
     334             :       END IF
     335     1004586 :       IF (PRESENT(drhob)) THEN
     336      572028 :          DO i = 1, 3
     337      429021 :             drhob(i)%array => rho_set%drhob(i)%array
     338      572028 :             CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_set%drhob(i)%array))
     339             :          END DO
     340             :       END IF
     341     1004586 :       IF (PRESENT(laplace_rhoa)) THEN
     342        3386 :          laplace_rhoa => rho_set%laplace_rhoa
     343        3386 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhoa))
     344             :       END IF
     345     1004586 :       IF (PRESENT(laplace_rhob)) THEN
     346        3386 :          laplace_rhob => rho_set%laplace_rhob
     347        3386 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(laplace_rhob))
     348             :       END IF
     349     1004586 :       IF (PRESENT(norm_drhoa)) THEN
     350      192845 :          norm_drhoa => rho_set%norm_drhoa
     351      192845 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhoa))
     352             :       END IF
     353     1004586 :       IF (PRESENT(norm_drhob)) THEN
     354      192845 :          norm_drhob => rho_set%norm_drhob
     355      192845 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(norm_drhob))
     356             :       END IF
     357     1004586 :       IF (PRESENT(rho_1_3)) THEN
     358       21173 :          rho_1_3 => rho_set%rho_1_3
     359       21173 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(rho_1_3))
     360             :       END IF
     361     1004586 :       IF (PRESENT(rhoa_1_3)) THEN
     362        3260 :          rhoa_1_3 => rho_set%rhoa_1_3
     363        3260 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(rhoa_1_3))
     364             :       END IF
     365     1004586 :       IF (PRESENT(rhob_1_3)) THEN
     366        3260 :          rhob_1_3 => rho_set%rhob_1_3
     367        3260 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(rhob_1_3))
     368             :       END IF
     369     1004586 :       IF (PRESENT(tau)) THEN
     370       14148 :          tau => rho_set%tau
     371       14148 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(tau))
     372             :       END IF
     373     1004586 :       IF (PRESENT(tau_a)) THEN
     374        3750 :          tau_a => rho_set%tau_a
     375        3750 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_a))
     376             :       END IF
     377     1004586 :       IF (PRESENT(tau_b)) THEN
     378        3750 :          tau_b => rho_set%tau_b
     379        3750 :          CPASSERT(my_can_return_null .OR. ASSOCIATED(tau_b))
     380             :       END IF
     381     1004586 :       IF (PRESENT(rho_cutoff)) rho_cutoff = rho_set%rho_cutoff
     382     1004586 :       IF (PRESENT(drho_cutoff)) drho_cutoff = rho_set%drho_cutoff
     383     1004586 :       IF (PRESENT(tau_cutoff)) tau_cutoff = rho_set%tau_cutoff
     384     2550526 :       IF (PRESENT(local_bounds)) local_bounds = rho_set%local_bounds
     385     1004586 :    END SUBROUTINE xc_rho_set_get
     386             : 
     387             :    #:mute
     388             :       #:def recover(variable)
     389             :          #! Determine component of the actual data
     390             :          #:set var_long_pw = (variable+"(i)" if variable.startswith("drho") else variable)
     391             :          #:set var_long_rho = (variable+"(i)%array" if variable.startswith("drho") else variable)
     392             :          #! Determine the flag name
     393             :          #! Remove spin states and potential underscore
     394             :          #:set is_1_3 = variable.endswith("_1_3")
     395             :          #:set var_base = variable.strip("_13")
     396             :          #:set is_spin = var_base.endswith("a") or var_base.endswith("b")
     397             :          #:set var_base = var_base.strip("ab_")
     398             :          #:set var_cflags = (var_base if not is_spin else var_base+"_spin")
     399             :          #:set var_cflags = (var_cflags if not is_1_3 else var_cflags+"_1_3")
     400             :          IF (PRESENT(${variable}$)) THEN
     401             :             #:if variable.startswith("drho")
     402             :             DO i = 1, 3
     403             :                #:else
     404             :                NULLIFY (${var_long_pw}$)
     405             :                ALLOCATE (${var_long_pw}$)
     406             :                #:endif
     407             :                CALL xc_rho_set_recover_pw_low(${var_long_pw}$, rho_set%${var_long_rho}$, pw_grid, pw_pool#{if variable =="drho"}#, rho_set%drhoa(i)%array, rho_set%drhob(i)%array#{endif}#)
     408             :                #:if not variable.startswith("drho")
     409             :                   NULLIFY (rho_set%${var_long_rho}$)
     410             :                #:else
     411             :                   END DO
     412             :                #:endif
     413             :                owns_data = #{if variable =="drho"}#.TRUE.#{else}#rho_set%owns%${var_cflags}$#{endif}#
     414             :             END IF
     415             :             #:enddef
     416             :          #:endmute
     417             : 
     418             : ! **************************************************************************************************
     419             : !> \brief Shifts association of the requested array to a pw grid
     420             : !>        Requires that the corresponding component of rho_set is associated
     421             : !>        If owns_data returns TRUE, the caller has to allocate the data later
     422             : !>        It is allowed to task for only one component per call.
     423             : !>        In case of drho, the array is allocated if not internally available and calculated from drhoa and drhob.
     424             : !> \param rho_set the object you want info about
     425             : !> \param pw_grid ...
     426             : !> \param pw_pool ...
     427             : !> \param owns_data ...
     428             : !> \param rho ...
     429             : !> \param drho ...
     430             : !> \param norm_drho ...
     431             : !> \param rhoa ...
     432             : !> \param rhob ...
     433             : !> \param norm_drhoa ...
     434             : !> \param norm_drhob ...
     435             : !> \param rho_1_3 ...
     436             : !> \param rhoa_1_3 ...
     437             : !> \param rhob_1_3 ...
     438             : !> \param laplace_rho ...
     439             : !> \param laplace_rhoa ...
     440             : !> \param laplace_rhob ...
     441             : !> \param drhoa ...
     442             : !> \param drhob ...
     443             : !> \param tau ...
     444             : !> \param tau_a ...
     445             : !> \param tau_b ...
     446             : ! **************************************************************************************************
     447      499038 :          SUBROUTINE xc_rho_set_recover_pw(rho_set, pw_grid, pw_pool, owns_data, rho, drho, norm_drho, &
     448             :                                           rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, &
     449             :                                           rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, &
     450             :                                           tau, tau_a, tau_b)
     451             :             TYPE(xc_rho_set_type)                              :: rho_set
     452             :             TYPE(pw_r3d_rs_type), DIMENSION(3), OPTIONAL, INTENT(OUT) :: drho, drhoa, drhob
     453             :             TYPE(pw_r3d_rs_type), OPTIONAL, POINTER                   :: rho, norm_drho, rhoa, rhob, norm_drhoa, &
     454             :                                                                          norm_drhob, rho_1_3, rhoa_1_3, &
     455             :                                                                          rhob_1_3, laplace_rho, laplace_rhoa, &
     456             :                                                                          laplace_rhob, tau, tau_a, tau_b
     457             :             TYPE(pw_grid_type), POINTER, INTENT(IN)            :: pw_grid
     458             :             TYPE(pw_pool_type), POINTER, INTENT(IN)            :: pw_pool
     459             :             LOGICAL, INTENT(OUT) :: owns_data
     460             : 
     461             :             INTEGER                                            :: i
     462             : 
     463             :             #:for variable in ["rho", "drho", "norm_drho", "rhoa", "rhob", "norm_drhoa", "norm_drhob", "rho_1_3", "rhoa_1_3", "rhob_1_3", "laplace_rho", "laplace_rhoa", "laplace_rhob", "drhoa", "drhob", "tau", "tau_a", "tau_b"]
     464      427405 :                $:recover(variable)
     465             :             #:endfor
     466             : 
     467       85481 :          END SUBROUTINE
     468             : 
     469      256443 :          SUBROUTINE xc_rho_set_recover_pw_low(rho_pw, rho, pw_grid, pw_pool, rhoa, rhob)
     470             :             TYPE(pw_r3d_rs_type), INTENT(OUT) :: rho_pw
     471             :             REAL(KIND=dp), DIMENSION(:, :, :), POINTER, CONTIGUOUS :: rho
     472             :             TYPE(pw_grid_type), POINTER, INTENT(IN) :: pw_grid
     473             :             TYPE(pw_pool_type), POINTER, INTENT(IN) :: pw_pool
     474             :             REAL(KIND=dp), DIMENSION(:, :, :), POINTER, OPTIONAL :: rhoa, rhob
     475             : 
     476      256443 :             IF (ASSOCIATED(rho)) THEN
     477      216981 :                CALL rho_pw%create(pw_grid=pw_grid, array_ptr=rho)
     478      216981 :                NULLIFY (rho)
     479       39462 :             ELSE IF (PRESENT(rhoa) .AND. PRESENT(rhob)) THEN
     480       39462 :                CALL pw_pool%create_pw(rho_pw)
     481       39462 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) SHARED(rho_pw,rhoa,rhob)
     482             :                rho_pw%array(:, :, :) = rhoa(:, :, :) + rhob(:, :, :)
     483             : !$OMP END PARALLEL WORKSHARE
     484             :             ELSE
     485             :                CALL cp_abort(__LOCATION__, "Either component or its spin parts (if applicable) "// &
     486           0 :                              "have to be associated in rho_set!")
     487             :             END IF
     488             : 
     489      256443 :          END SUBROUTINE
     490             : 
     491             : ! **************************************************************************************************
     492             : !> \brief cleans (releases) most of the data stored in the rho_set giving back
     493             : !>      what it can to the pw_pool
     494             : !> \param rho_set the rho_set to be cleaned
     495             : !> \param pw_pool place to give back 3d arrays,...
     496             : !> \author Fawzi Mohamed
     497             : ! **************************************************************************************************
     498      269750 :          SUBROUTINE xc_rho_set_clean(rho_set, pw_pool)
     499             :             TYPE(xc_rho_set_type), INTENT(INOUT)               :: rho_set
     500             :             TYPE(pw_pool_type), POINTER                        :: pw_pool
     501             : 
     502             :             INTEGER                                            :: idir
     503             : 
     504      269750 :             IF (rho_set%owns%rho) THEN
     505      241864 :                CALL pw_pool%give_back_cr3d(rho_set%rho)
     506             :             ELSE
     507       27886 :                NULLIFY (rho_set%rho)
     508             :             END IF
     509      269750 :             IF (rho_set%owns%rho_1_3) THEN
     510      151871 :                CALL pw_pool%give_back_cr3d(rho_set%rho_1_3)
     511             :             ELSE
     512      117879 :                NULLIFY (rho_set%rho_1_3)
     513             :             END IF
     514      269750 :             IF (rho_set%owns%drho) THEN
     515      772936 :                DO idir = 1, 3
     516      772936 :                   CALL pw_pool%give_back_cr3d(rho_set%drho(idir)%array)
     517             :                END DO
     518             :             ELSE
     519      306064 :                DO idir = 1, 3
     520      306064 :                   NULLIFY (rho_set%drho(idir)%array)
     521             :                END DO
     522             :             END IF
     523      269750 :             IF (rho_set%owns%norm_drho) THEN
     524      211840 :                CALL pw_pool%give_back_cr3d(rho_set%norm_drho)
     525             :             ELSE
     526       57910 :                NULLIFY (rho_set%norm_drho)
     527             :             END IF
     528      269750 :             IF (rho_set%owns%laplace_rho) THEN
     529      143761 :                CALL pw_pool%give_back_cr3d(rho_set%laplace_rho)
     530             :             ELSE
     531      125989 :                NULLIFY (rho_set%laplace_rho)
     532             :             END IF
     533      269750 :             IF (rho_set%owns%tau) THEN
     534      143111 :                CALL pw_pool%give_back_cr3d(rho_set%tau)
     535             :             ELSE
     536      126639 :                NULLIFY (rho_set%tau)
     537             :             END IF
     538      269750 :             IF (rho_set%owns%rho_spin) THEN
     539      170997 :                CALL pw_pool%give_back_cr3d(rho_set%rhoa)
     540      170997 :                CALL pw_pool%give_back_cr3d(rho_set%rhob)
     541             :             ELSE
     542       98753 :                NULLIFY (rho_set%rhoa, rho_set%rhob)
     543             :             END IF
     544      269750 :             IF (rho_set%owns%rho_spin_1_3) THEN
     545      144729 :                CALL pw_pool%give_back_cr3d(rho_set%rhoa_1_3)
     546      144729 :                CALL pw_pool%give_back_cr3d(rho_set%rhob_1_3)
     547             :             ELSE
     548      125021 :                NULLIFY (rho_set%rhoa_1_3, rho_set%rhob_1_3)
     549             :             END IF
     550      269750 :             IF (rho_set%owns%drho_spin) THEN
     551      631148 :                DO idir = 1, 3
     552      473361 :                   CALL pw_pool%give_back_cr3d(rho_set%drhoa(idir)%array)
     553      631148 :                   CALL pw_pool%give_back_cr3d(rho_set%drhob(idir)%array)
     554             :                END DO
     555             :             ELSE
     556      447852 :                DO idir = 1, 3
     557      447852 :                   NULLIFY (rho_set%drhoa(idir)%array, rho_set%drhob(idir)%array)
     558             :                END DO
     559             :             END IF
     560      269750 :             IF (rho_set%owns%laplace_rho_spin) THEN
     561      143321 :                CALL pw_pool%give_back_cr3d(rho_set%laplace_rhoa)
     562      143321 :                CALL pw_pool%give_back_cr3d(rho_set%laplace_rhob)
     563             :             ELSE
     564      126429 :                NULLIFY (rho_set%laplace_rhoa, rho_set%laplace_rhob)
     565             :             END IF
     566      269750 :             IF (rho_set%owns%norm_drho_spin) THEN
     567      160133 :                CALL pw_pool%give_back_cr3d(rho_set%norm_drhoa)
     568      160133 :                CALL pw_pool%give_back_cr3d(rho_set%norm_drhob)
     569             :             ELSE
     570      109617 :                NULLIFY (rho_set%norm_drhoa, rho_set%norm_drhob)
     571             :             END IF
     572      269750 :             IF (rho_set%owns%tau_spin) THEN
     573      143111 :                CALL pw_pool%give_back_cr3d(rho_set%tau_a)
     574      143111 :                CALL pw_pool%give_back_cr3d(rho_set%tau_b)
     575             :             ELSE
     576      126639 :                NULLIFY (rho_set%tau_a, rho_set%tau_b)
     577             :             END IF
     578             : 
     579      269750 :             CALL xc_rho_cflags_setall(rho_set%has, .FALSE.)
     580      269750 :             CALL xc_rho_cflags_setall(rho_set%owns, .FALSE.)
     581             : 
     582      269750 :          END SUBROUTINE xc_rho_set_clean
     583             : 
     584             : ! **************************************************************************************************
     585             : !> \brief updates the given rho set with the density given by
     586             : !>      rho_r (and rho_g). The rho set will contain the components specified
     587             : !>      in needs
     588             : !> \param rho_set the rho_set to update
     589             : !> \param rho_r the new density (in r space)
     590             : !> \param rho_g the new density (in g space, needed for some
     591             : !>        derivatives)
     592             : !> \param tau ...
     593             : !> \param needs the components of rho that are needed
     594             : !> \param xc_deriv_method_id ...
     595             : !> \param xc_rho_smooth_id ...
     596             : !> \param pw_pool pool for the allocation of pw and cr3d
     597             : ! **************************************************************************************************
     598      150683 :          SUBROUTINE xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
     599             :                                       xc_deriv_method_id, xc_rho_smooth_id, pw_pool)
     600             :             TYPE(xc_rho_set_type), INTENT(INOUT)               :: rho_set
     601             :             TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN)            :: rho_r
     602             :             TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER               :: rho_g
     603             :             TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER, INTENT(IN)   :: tau
     604             :             TYPE(xc_rho_cflags_type), INTENT(in)               :: needs
     605             :             INTEGER, INTENT(IN)                                :: xc_deriv_method_id, xc_rho_smooth_id
     606             :             TYPE(pw_pool_type), POINTER                        :: pw_pool
     607             : 
     608             :             REAL(KIND=dp), PARAMETER                           :: f13 = (1.0_dp/3.0_dp)
     609             : 
     610             :             INTEGER                                            :: i, idir, ispin, j, k, nspins
     611             :             LOGICAL                                            :: gradient_f, my_rho_g_local, &
     612             :                                                                   needs_laplace, needs_rho_g
     613             :             REAL(kind=dp)                                      :: rho_cutoff
     614      452049 :             TYPE(pw_r3d_rs_type), DIMENSION(2)                      :: laplace_rho_r
     615     1356147 :             TYPE(pw_r3d_rs_type), DIMENSION(3, 2)                   :: drho_r
     616             :             TYPE(pw_c1d_gs_type)                                      :: my_rho_g, tmp_g
     617      452049 :             TYPE(pw_r3d_rs_type), DIMENSION(2)                        :: my_rho_r
     618             : 
     619     1506830 :             IF (ANY(rho_set%local_bounds /= pw_pool%pw_grid%bounds_local)) &
     620           0 :                CPABORT("pw_pool cr3d have different size than expected")
     621      150683 :             nspins = SIZE(rho_r)
     622     1506830 :             rho_set%local_bounds = rho_r(1)%pw_grid%bounds_local
     623      150683 :             rho_cutoff = 0.5*rho_set%rho_cutoff
     624             : 
     625      150683 :             my_rho_g_local = .FALSE.
     626             :             ! some checks
     627      119655 :             SELECT CASE (nspins)
     628             :             CASE (1)
     629      119655 :                CPASSERT(.NOT. needs%rho_spin)
     630      119655 :                CPASSERT(.NOT. needs%drho_spin)
     631      119655 :                CPASSERT(.NOT. needs%norm_drho_spin)
     632      119655 :                CPASSERT(.NOT. needs%rho_spin_1_3)
     633      119655 :                CPASSERT(.NOT. needs%tau_spin)
     634      119655 :                CPASSERT(.NOT. needs%laplace_rho_spin)
     635             :             CASE (2)
     636       31028 :                CPASSERT(.NOT. needs%rho)
     637       31028 :                CPASSERT(.NOT. needs%drho)
     638       31028 :                CPASSERT(.NOT. needs%rho_1_3)
     639       31028 :                CPASSERT(.NOT. needs%tau)
     640       31028 :                CPASSERT(.NOT. needs%laplace_rho)
     641             :             CASE default
     642      150683 :                CPABORT("Unknown number of spin states")
     643             :             END SELECT
     644             : 
     645      150683 :             CALL xc_rho_set_clean(rho_set, pw_pool=pw_pool)
     646             : 
     647      150683 :             needs_laplace = (needs%laplace_rho .OR. needs%laplace_rho_spin)
     648             :             gradient_f = (needs%drho_spin .OR. needs%norm_drho_spin .OR. &
     649             :                           needs%drho .OR. needs%norm_drho .OR. &
     650      150683 :                           needs_laplace)
     651             :             needs_rho_g = ((xc_deriv_method_id == xc_deriv_spline3 .OR. &
     652             :                             xc_deriv_method_id == xc_deriv_spline2 .OR. &
     653      150683 :                             xc_deriv_method_id == xc_deriv_pw)) .AND. (gradient_f .OR. needs_laplace)
     654      150683 :             IF ((gradient_f .AND. needs_laplace) .AND. &
     655             :                 (xc_deriv_method_id /= xc_deriv_pw)) THEN
     656             :                CALL cp_abort(__LOCATION__, &
     657             :                              "MGGA functionals that require the Laplacian are "// &
     658           0 :                              "only compatible with 'XC_DERIV PW' and 'XC_SMOOTH_RHO NONE'")
     659             :             END IF
     660             : 
     661      150683 :             IF (needs_rho_g) THEN
     662       81409 :                CALL pw_pool%create_pw(tmp_g)
     663             :             END IF
     664      332394 :             DO ispin = 1, nspins
     665      181711 :                CALL pw_pool%create_pw(my_rho_r(ispin))
     666             :                ! introduce a smoothing kernel on the density
     667      181711 :                IF (xc_rho_smooth_id == xc_rho_no_smooth) THEN
     668      181213 :                   IF (needs_rho_g) THEN
     669       99231 :                      IF (ASSOCIATED(rho_g)) THEN
     670       83899 :                         my_rho_g_local = .FALSE.
     671       83899 :                         my_rho_g = rho_g(ispin)
     672             :                      END IF
     673             :                   END IF
     674             : 
     675      181213 :                   CALL pw_copy(rho_r(ispin), my_rho_r(ispin))
     676             :                ELSE
     677         498 :                   CALL xc_pw_smooth(rho_r(ispin), my_rho_r(ispin), xc_rho_smooth_id)
     678             :                END IF
     679             : 
     680      332394 :                IF (gradient_f) THEN ! calculate the grad of rho
     681             :                   ! normally when you need the gradient you need the whole gradient
     682             :                   ! (for the partial integration)
     683             :                   ! deriv rho
     684      408004 :                   DO idir = 1, 3
     685      408004 :                      CALL pw_pool%create_pw(drho_r(idir, ispin))
     686             :                   END DO
     687      102001 :                   IF (needs_rho_g) THEN
     688       99297 :                      IF (.NOT. ASSOCIATED(my_rho_g%pw_grid)) THEN
     689       15398 :                         my_rho_g_local = .TRUE.
     690       15398 :                         CALL pw_pool%create_pw(my_rho_g)
     691       15398 :                         CALL pw_transfer(my_rho_r(ispin), my_rho_g)
     692             :                      END IF
     693       83899 :                      IF (.NOT. my_rho_g_local .AND. (xc_deriv_method_id == xc_deriv_spline2 .OR. &
     694             :                                                      xc_deriv_method_id == xc_deriv_spline3)) THEN
     695        7140 :                         CALL pw_pool%create_pw(my_rho_g)
     696        7140 :                         my_rho_g_local = .TRUE.
     697        7140 :                         CALL pw_copy(rho_g(ispin), my_rho_g)
     698             :                      END IF
     699             :                   END IF
     700      102001 :                   IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
     701        1332 :                      CALL pw_pool%create_pw(laplace_rho_r(ispin))
     702        1332 :                      CALL xc_pw_laplace(my_rho_g, pw_pool, xc_deriv_method_id, laplace_rho_r(ispin), tmp_g=tmp_g)
     703             :                   END IF
     704      102001 :                   CALL xc_pw_gradient(my_rho_r(ispin), my_rho_g, tmp_g, drho_r(:, ispin), xc_deriv_method_id)
     705             : 
     706      102001 :                   IF (needs_rho_g) THEN
     707       99297 :                      IF (my_rho_g_local) THEN
     708       22538 :                         my_rho_g_local = .FALSE.
     709       22538 :                         CALL pw_pool%give_back_pw(my_rho_g)
     710             :                      END IF
     711             :                   END IF
     712             : 
     713      102001 :                   IF (xc_deriv_method_id /= xc_deriv_pw) THEN
     714        9898 :                      CALL pw_spline_scale_deriv(drho_r(:, ispin))
     715             :                   END IF
     716             : 
     717             :                END IF
     718             : 
     719             :             END DO
     720             : 
     721      150683 :             IF (ASSOCIATED(tmp_g%pw_grid)) THEN
     722       81409 :                CALL pw_pool%give_back_pw(tmp_g)
     723             :             END IF
     724             : 
     725      119655 :             SELECT CASE (nspins)
     726             :             CASE (1)
     727      119655 :                IF (needs%rho_1_3) THEN
     728        9986 :                   CALL pw_pool%create_cr3d(rho_set%rho_1_3)
     729        9986 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
     730             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     731             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     732             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     733             :                            rho_set%rho_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
     734             :                         END DO
     735             :                      END DO
     736             :                   END DO
     737        9986 :                   rho_set%owns%rho_1_3 = .TRUE.
     738        9986 :                   rho_set%has%rho_1_3 = .TRUE.
     739             :                END IF
     740      119655 :                IF (needs%rho) THEN
     741      119655 :                   rho_set%rho => my_rho_r(1)%array
     742      119655 :                   NULLIFY (my_rho_r(1)%array)
     743      119655 :                   rho_set%owns%rho = .TRUE.
     744      119655 :                   rho_set%has%rho = .TRUE.
     745             :                END IF
     746      119655 :                IF (needs%norm_drho) THEN
     747       64925 :                   CALL pw_pool%create_cr3d(rho_set%norm_drho)
     748       64925 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
     749             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     750             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     751             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     752             :                            rho_set%norm_drho(i, j, k) = SQRT( &
     753             :                                                         drho_r(1, 1)%array(i, j, k)**2 + &
     754             :                                                         drho_r(2, 1)%array(i, j, k)**2 + &
     755             :                                                         drho_r(3, 1)%array(i, j, k)**2)
     756             :                         END DO
     757             :                      END DO
     758             :                   END DO
     759       64925 :                   rho_set%owns%norm_drho = .TRUE.
     760       64925 :                   rho_set%has%norm_drho = .TRUE.
     761             :                END IF
     762      119655 :                IF (needs%laplace_rho) THEN
     763         840 :                   rho_set%laplace_rho => laplace_rho_r(1)%array
     764         840 :                   NULLIFY (laplace_rho_r(1)%array)
     765         840 :                   rho_set%owns%laplace_rho = .TRUE.
     766         840 :                   rho_set%has%laplace_rho = .TRUE.
     767             :                END IF
     768             : 
     769      119655 :                IF (needs%drho) THEN
     770      249196 :                   DO idir = 1, 3
     771      186897 :                      rho_set%drho(idir)%array => drho_r(idir, 1)%array
     772      249196 :                      NULLIFY (drho_r(idir, 1)%array)
     773             :                   END DO
     774       62299 :                   rho_set%owns%drho = .TRUE.
     775       62299 :                   rho_set%has%drho = .TRUE.
     776             :                END IF
     777             :             CASE (2)
     778       31028 :                IF (needs%rho_spin_1_3) THEN
     779        1630 :                   CALL pw_pool%create_cr3d(rho_set%rhoa_1_3)
     780             :                   !assume that the bounds are the same?
     781        1630 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
     782             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     783             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     784             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     785             :                            rho_set%rhoa_1_3(i, j, k) = MAX(my_rho_r(1)%array(i, j, k), 0.0_dp)**f13
     786             :                         END DO
     787             :                      END DO
     788             :                   END DO
     789        1630 :                   CALL pw_pool%create_cr3d(rho_set%rhob_1_3)
     790             :                   !assume that the bounds are the same?
     791        1630 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,my_rho_r)
     792             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     793             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     794             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     795             :                            rho_set%rhob_1_3(i, j, k) = MAX(my_rho_r(2)%array(i, j, k), 0.0_dp)**f13
     796             :                         END DO
     797             :                      END DO
     798             :                   END DO
     799        1630 :                   rho_set%owns%rho_spin_1_3 = .TRUE.
     800        1630 :                   rho_set%has%rho_spin_1_3 = .TRUE.
     801             :                END IF
     802       31028 :                IF (needs%norm_drho) THEN
     803             : 
     804       17476 :                   CALL pw_pool%create_cr3d(rho_set%norm_drho)
     805       17476 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
     806             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     807             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     808             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     809             :                            rho_set%norm_drho(i, j, k) = SQRT( &
     810             :                                                         (drho_r(1, 1)%array(i, j, k) + drho_r(1, 2)%array(i, j, k))**2 + &
     811             :                                                         (drho_r(2, 1)%array(i, j, k) + drho_r(2, 2)%array(i, j, k))**2 + &
     812             :                                                         (drho_r(3, 1)%array(i, j, k) + drho_r(3, 2)%array(i, j, k))**2)
     813             :                         END DO
     814             :                      END DO
     815             :                   END DO
     816             : 
     817       17476 :                   rho_set%owns%norm_drho = .TRUE.
     818       17476 :                   rho_set%has%norm_drho = .TRUE.
     819             :                END IF
     820       31028 :                IF (needs%rho_spin) THEN
     821             : 
     822       31028 :                   rho_set%rhoa => my_rho_r(1)%array
     823       31028 :                   NULLIFY (my_rho_r(1)%array)
     824             : 
     825       31028 :                   rho_set%rhob => my_rho_r(2)%array
     826       31028 :                   NULLIFY (my_rho_r(2)%array)
     827             : 
     828       31028 :                   rho_set%owns%rho_spin = .TRUE.
     829       31028 :                   rho_set%has%rho_spin = .TRUE.
     830             :                END IF
     831       31028 :                IF (needs%norm_drho_spin) THEN
     832             : 
     833       18538 :                   CALL pw_pool%create_cr3d(rho_set%norm_drhoa)
     834       18538 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
     835             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     836             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     837             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     838             :                            rho_set%norm_drhoa(i, j, k) = SQRT( &
     839             :                                                          drho_r(1, 1)%array(i, j, k)**2 + &
     840             :                                                          drho_r(2, 1)%array(i, j, k)**2 + &
     841             :                                                          drho_r(3, 1)%array(i, j, k)**2)
     842             :                         END DO
     843             :                      END DO
     844             :                   END DO
     845             : 
     846       18538 :                   CALL pw_pool%create_cr3d(rho_set%norm_drhob)
     847       18538 :                   rho_set%owns%norm_drho_spin = .TRUE.
     848       18538 : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k) SHARED(rho_set,drho_r)
     849             :                   DO k = rho_set%local_bounds(1, 3), rho_set%local_bounds(2, 3)
     850             :                      DO j = rho_set%local_bounds(1, 2), rho_set%local_bounds(2, 2)
     851             :                         DO i = rho_set%local_bounds(1, 1), rho_set%local_bounds(2, 1)
     852             :                            rho_set%norm_drhob(i, j, k) = SQRT( &
     853             :                                                          drho_r(1, 2)%array(i, j, k)**2 + &
     854             :                                                          drho_r(2, 2)%array(i, j, k)**2 + &
     855             :                                                          drho_r(3, 2)%array(i, j, k)**2)
     856             :                         END DO
     857             :                      END DO
     858             :                   END DO
     859             : 
     860       18538 :                   rho_set%owns%norm_drho_spin = .TRUE.
     861       18538 :                   rho_set%has%norm_drho_spin = .TRUE.
     862             :                END IF
     863       31028 :                IF (needs%laplace_rho_spin) THEN
     864         246 :                   rho_set%laplace_rhoa => laplace_rho_r(1)%array
     865         246 :                   NULLIFY (laplace_rho_r(1)%array)
     866             : 
     867         246 :                   rho_set%laplace_rhob => laplace_rho_r(2)%array
     868         246 :                   NULLIFY (laplace_rho_r(2)%array)
     869             : 
     870         246 :                   rho_set%owns%laplace_rho_spin = .TRUE.
     871         246 :                   rho_set%has%laplace_rho_spin = .TRUE.
     872             :                END IF
     873      181711 :                IF (needs%drho_spin) THEN
     874       64768 :                   DO idir = 1, 3
     875       48576 :                      rho_set%drhoa(idir)%array => drho_r(idir, 1)%array
     876       48576 :                      NULLIFY (drho_r(idir, 1)%array)
     877       48576 :                      rho_set%drhob(idir)%array => drho_r(idir, 2)%array
     878       64768 :                      NULLIFY (drho_r(idir, 2)%array)
     879             :                   END DO
     880       16192 :                   rho_set%owns%drho_spin = .TRUE.
     881       16192 :                   rho_set%has%drho_spin = .TRUE.
     882             :                END IF
     883             :             END SELECT
     884             :             ! post cleanup
     885      332394 :             DO ispin = 1, nspins
     886      181711 :                IF (needs%laplace_rho .OR. needs%laplace_rho_spin) THEN
     887        1332 :                   CALL pw_pool%give_back_pw(laplace_rho_r(ispin))
     888             :                END IF
     889      877527 :                DO idir = 1, 3
     890      726844 :                   CALL pw_pool%give_back_pw(drho_r(idir, ispin))
     891             :                END DO
     892             :             END DO
     893      332394 :             DO ispin = 1, nspins
     894      332394 :                CALL pw_pool%give_back_pw(my_rho_r(ispin))
     895             :             END DO
     896             : 
     897             :             ! tau part
     898      150683 :             IF (needs%tau .OR. needs%tau_spin) THEN
     899        3480 :                CPASSERT(ASSOCIATED(tau))
     900        7668 :                DO ispin = 1, nspins
     901      154871 :                   CPASSERT(ASSOCIATED(tau(ispin)%array))
     902             :                END DO
     903             :             END IF
     904      150683 :             IF (needs%tau) THEN
     905        2772 :                rho_set%tau => tau(1)%array
     906        2772 :                rho_set%owns%tau = .FALSE.
     907        2772 :                rho_set%has%tau = .TRUE.
     908             :             END IF
     909      150683 :             IF (needs%tau_spin) THEN
     910         708 :                rho_set%tau_a => tau(1)%array
     911         708 :                rho_set%tau_b => tau(2)%array
     912         708 :                rho_set%owns%tau_spin = .FALSE.
     913         708 :                rho_set%has%tau_spin = .TRUE.
     914             :             END IF
     915             : 
     916      150683 :             CPASSERT(xc_rho_cflags_equal(rho_set%has, needs))
     917             : 
     918      301366 :          END SUBROUTINE xc_rho_set_update
     919             : 
     920           0 :       END MODULE xc_rho_set_types

Generated by: LCOV version 1.15