LCOV - code coverage report
Current view: top level - src/xc - xc_tfw.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 0.0 % 153 0
Test Date: 2025-07-25 12:55:17 Functions: 0.0 % 14 0

            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 Calculate the Thomas-Fermi kinetic energy functional
      10              : !>      plus the von Weizsaecker term
      11              : !> \par History
      12              : !>      JGH (26.02.2003) : OpenMP enabled
      13              : !>      fawzi (04.2004)  : adapted to the new xc interface
      14              : !> \author JGH (18.02.2002)
      15              : ! **************************************************************************************************
      16              : MODULE xc_tfw
      17              :    USE cp_array_utils,                  ONLY: cp_3d_r_cp_type
      18              :    USE kinds,                           ONLY: dp
      19              :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      20              :                                               deriv_norm_drhoa,&
      21              :                                               deriv_norm_drhob,&
      22              :                                               deriv_rho,&
      23              :                                               deriv_rhoa,&
      24              :                                               deriv_rhob
      25              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      26              :                                               xc_dset_get_derivative
      27              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      28              :                                               xc_derivative_type
      29              :    USE xc_functionals_utilities,        ONLY: set_util
      30              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      31              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      32              :                                               xc_rho_set_type
      33              : #include "../base/base_uses.f90"
      34              : 
      35              :    IMPLICIT NONE
      36              : 
      37              :    PRIVATE
      38              : 
      39              : ! *** Global parameters ***
      40              : 
      41              :    REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
      42              :    REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
      43              :                                f23 = 2.0_dp*f13, &
      44              :                                f43 = 4.0_dp*f13, &
      45              :                                f53 = 5.0_dp*f13
      46              : 
      47              :    PUBLIC :: tfw_lda_info, tfw_lda_eval, tfw_lsd_info, tfw_lsd_eval
      48              : 
      49              :    REAL(KIND=dp) :: cf, flda, flsd, fvw
      50              :    REAL(KIND=dp) :: eps_rho
      51              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tfw'
      52              : 
      53              : CONTAINS
      54              : 
      55              : ! **************************************************************************************************
      56              : !> \brief ...
      57              : !> \param cutoff ...
      58              : ! **************************************************************************************************
      59            0 :    SUBROUTINE tfw_init(cutoff)
      60              : 
      61              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      62              : 
      63            0 :       eps_rho = cutoff
      64            0 :       CALL set_util(cutoff)
      65              : 
      66            0 :       cf = 0.3_dp*(3.0_dp*pi*pi)**f23
      67            0 :       flda = cf
      68            0 :       flsd = flda*2.0_dp**f23
      69            0 :       fvw = 1.0_dp/72.0_dp
      70              : 
      71            0 :    END SUBROUTINE tfw_init
      72              : 
      73              : ! **************************************************************************************************
      74              : !> \brief ...
      75              : !> \param reference ...
      76              : !> \param shortform ...
      77              : !> \param needs ...
      78              : !> \param max_deriv ...
      79              : ! **************************************************************************************************
      80            0 :    SUBROUTINE tfw_lda_info(reference, shortform, needs, max_deriv)
      81              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      82              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      83              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      84              : 
      85            0 :       IF (PRESENT(reference)) THEN
      86            0 :          reference = "Thomas-Fermi-Weizsaecker kinetic energy functional {LDA version}"
      87              :       END IF
      88            0 :       IF (PRESENT(shortform)) THEN
      89            0 :          shortform = "TF+vW kinetic energy functional {LDA}"
      90              :       END IF
      91            0 :       IF (PRESENT(needs)) THEN
      92            0 :          needs%rho = .TRUE.
      93            0 :          needs%rho_1_3 = .TRUE.
      94            0 :          needs%norm_drho = .TRUE.
      95              :       END IF
      96            0 :       IF (PRESENT(max_deriv)) max_deriv = 3
      97              : 
      98            0 :    END SUBROUTINE tfw_lda_info
      99              : 
     100              : ! **************************************************************************************************
     101              : !> \brief ...
     102              : !> \param reference ...
     103              : !> \param shortform ...
     104              : !> \param needs ...
     105              : !> \param max_deriv ...
     106              : ! **************************************************************************************************
     107            0 :    SUBROUTINE tfw_lsd_info(reference, shortform, needs, max_deriv)
     108              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
     109              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
     110              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
     111              : 
     112            0 :       IF (PRESENT(reference)) THEN
     113            0 :          reference = "Thomas-Fermi-Weizsaecker kinetic energy functional"
     114              :       END IF
     115            0 :       IF (PRESENT(shortform)) THEN
     116            0 :          shortform = "TF+vW kinetic energy functional"
     117              :       END IF
     118            0 :       IF (PRESENT(needs)) THEN
     119            0 :          needs%rho_spin = .TRUE.
     120            0 :          needs%rho_spin_1_3 = .TRUE.
     121            0 :          needs%norm_drho = .TRUE.
     122              :       END IF
     123            0 :       IF (PRESENT(max_deriv)) max_deriv = 3
     124              : 
     125            0 :    END SUBROUTINE tfw_lsd_info
     126              : 
     127              : ! **************************************************************************************************
     128              : !> \brief ...
     129              : !> \param rho_set ...
     130              : !> \param deriv_set ...
     131              : !> \param order ...
     132              : ! **************************************************************************************************
     133            0 :    SUBROUTINE tfw_lda_eval(rho_set, deriv_set, order)
     134              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     135              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     136              :       INTEGER, INTENT(in)                                :: order
     137              : 
     138              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'tfw_lda_eval'
     139              : 
     140              :       INTEGER                                            :: handle, npoints
     141              :       INTEGER, DIMENSION(2, 3)                           :: bo
     142              :       REAL(KIND=dp)                                      :: epsilon_rho
     143            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: s
     144            0 :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
     145            0 :          e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, e_rho_rho_rho, grho, &
     146            0 :          r13, rho
     147              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     148              : 
     149            0 :       CALL timeset(routineN, handle)
     150              : 
     151              :       CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
     152            0 :                           norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho)
     153            0 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     154            0 :       CALL tfw_init(epsilon_rho)
     155              : 
     156            0 :       ALLOCATE (s(npoints))
     157            0 :       CALL calc_s(rho, grho, s, npoints)
     158              : 
     159            0 :       IF (order >= 0) THEN
     160              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     161            0 :                                          allocate_deriv=.TRUE.)
     162            0 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     163              : 
     164            0 :          CALL tfw_u_0(rho, r13, s, e_0, npoints)
     165              :       END IF
     166            0 :       IF (order >= 1 .OR. order == -1) THEN
     167              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     168            0 :                                          allocate_deriv=.TRUE.)
     169            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     170              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     171            0 :                                          allocate_deriv=.TRUE.)
     172            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     173              : 
     174            0 :          CALL tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
     175              :       END IF
     176            0 :       IF (order >= 2 .OR. order == -2) THEN
     177              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     178            0 :                                          allocate_deriv=.TRUE.)
     179            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     180              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
     181            0 :                                          allocate_deriv=.TRUE.)
     182            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
     183              :          deriv => xc_dset_get_derivative(deriv_set, &
     184            0 :                                          [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     185            0 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     186              : 
     187              :          CALL tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, &
     188            0 :                       e_ndrho_ndrho, npoints)
     189              :       END IF
     190            0 :       IF (order >= 3 .OR. order == -3) THEN
     191              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     192            0 :                                          allocate_deriv=.TRUE.)
     193            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     194              :          deriv => xc_dset_get_derivative(deriv_set, &
     195            0 :                                          [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
     196            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
     197              :          deriv => xc_dset_get_derivative(deriv_set, &
     198            0 :                                          [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
     199            0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
     200              : 
     201              :          CALL tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
     202            0 :                       e_rho_ndrho_ndrho, npoints)
     203              :       END IF
     204            0 :       IF (order > 3 .OR. order < -3) THEN
     205            0 :          CPABORT("derivatives bigger than 3 not implemented")
     206              :       END IF
     207              : 
     208            0 :       DEALLOCATE (s)
     209            0 :       CALL timestop(handle)
     210            0 :    END SUBROUTINE tfw_lda_eval
     211              : 
     212              : ! **************************************************************************************************
     213              : !> \brief ...
     214              : !> \param rho ...
     215              : !> \param grho ...
     216              : !> \param s ...
     217              : !> \param npoints ...
     218              : ! **************************************************************************************************
     219            0 :    SUBROUTINE calc_s(rho, grho, s, npoints)
     220              :       REAL(KIND=dp), DIMENSION(*), INTENT(in)            :: rho, grho
     221              :       REAL(KIND=dp), DIMENSION(*), INTENT(out)           :: s
     222              :       INTEGER, INTENT(in)                                :: npoints
     223              : 
     224              :       INTEGER                                            :: ip
     225              : 
     226              : !$OMP     PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     227            0 : !$OMP     SHARED(npoints,rho,eps_rho,s,grho)
     228              :       DO ip = 1, npoints
     229              :          IF (rho(ip) < eps_rho) THEN
     230              :             s(ip) = 0.0_dp
     231              :          ELSE
     232              :             s(ip) = grho(ip)*grho(ip)/rho(ip)
     233              :          END IF
     234              :       END DO
     235            0 :    END SUBROUTINE calc_s
     236              : 
     237              : ! **************************************************************************************************
     238              : !> \brief ...
     239              : !> \param rho_set ...
     240              : !> \param deriv_set ...
     241              : !> \param order ...
     242              : ! **************************************************************************************************
     243            0 :    SUBROUTINE tfw_lsd_eval(rho_set, deriv_set, order)
     244              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     245              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     246              :       INTEGER, INTENT(in)                                :: order
     247              : 
     248              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'tfw_lsd_eval'
     249              :       INTEGER, DIMENSION(2), PARAMETER :: &
     250              :          norm_drho_spin_name = [deriv_norm_drhoa, deriv_norm_drhob], &
     251              :          rho_spin_name = [deriv_rhoa, deriv_rhob]
     252              : 
     253              :       INTEGER                                            :: handle, i, ispin, npoints
     254              :       INTEGER, DIMENSION(2, 3)                           :: bo
     255              :       REAL(KIND=dp)                                      :: epsilon_rho
     256            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: s
     257              :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     258            0 :          POINTER                                         :: e_0, e_ndrho, e_ndrho_ndrho, e_rho, &
     259            0 :                                                             e_rho_ndrho, e_rho_ndrho_ndrho, &
     260            0 :                                                             e_rho_rho, e_rho_rho_ndrho, &
     261            0 :                                                             e_rho_rho_rho
     262            0 :       TYPE(cp_3d_r_cp_type), DIMENSION(2)                :: norm_drho, rho, rho_1_3
     263              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     264              : 
     265            0 :       CALL timeset(routineN, handle)
     266            0 :       NULLIFY (deriv)
     267            0 :       DO i = 1, 2
     268            0 :          NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
     269              :       END DO
     270              : 
     271              :       CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
     272              :                           rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
     273              :                           rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
     274              :                           norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
     275            0 :                           local_bounds=bo)
     276            0 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     277            0 :       CALL tfw_init(epsilon_rho)
     278              : 
     279            0 :       ALLOCATE (s(npoints))
     280              : 
     281            0 :       DO ispin = 1, 2
     282            0 :          CALL calc_s(rho(ispin)%array, norm_drho(ispin)%array, s, npoints)
     283              : 
     284            0 :          IF (order >= 0) THEN
     285              :             deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     286            0 :                                             allocate_deriv=.TRUE.)
     287            0 :             CALL xc_derivative_get(deriv, deriv_data=e_0)
     288              : 
     289              :             CALL tfw_p_0(rho(ispin)%array, &
     290            0 :                          rho_1_3(ispin)%array, s, e_0, npoints)
     291              :          END IF
     292            0 :          IF (order >= 1 .OR. order == -1) THEN
     293              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
     294            0 :                                             allocate_deriv=.TRUE.)
     295            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho)
     296              :             deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin)], &
     297            0 :                                             allocate_deriv=.TRUE.)
     298            0 :             CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     299              : 
     300              :             CALL tfw_p_1(rho(ispin)%array, norm_drho(ispin)%array, &
     301            0 :                          rho_1_3(ispin)%array, s, e_rho, e_ndrho, npoints)
     302              :          END IF
     303            0 :          IF (order >= 2 .OR. order == -2) THEN
     304              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     305            0 :                                                         rho_spin_name(ispin)], allocate_deriv=.TRUE.)
     306            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     307              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     308            0 :                                                         norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
     309            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
     310              :             deriv => xc_dset_get_derivative(deriv_set, [norm_drho_spin_name(ispin), &
     311            0 :                                                         norm_drho_spin_name(ispin)], allocate_deriv=.TRUE.)
     312            0 :             CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
     313              : 
     314              :             CALL tfw_p_2(rho(ispin)%array, norm_drho(ispin)%array, &
     315              :                          rho_1_3(ispin)%array, s, e_rho_rho, e_rho_ndrho, &
     316            0 :                          e_ndrho_ndrho, npoints)
     317              :          END IF
     318            0 :          IF (order >= 3 .OR. order == -3) THEN
     319              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     320              :                                                         rho_spin_name(ispin), rho_spin_name(ispin)], &
     321            0 :                                             allocate_deriv=.TRUE.)
     322            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     323              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     324              :                                                         rho_spin_name(ispin), norm_drho_spin_name(ispin)], &
     325            0 :                                             allocate_deriv=.TRUE.)
     326            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
     327              :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     328              :                                                         norm_drho_spin_name(ispin), norm_drho_spin_name(ispin)], &
     329            0 :                                             allocate_deriv=.TRUE.)
     330            0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
     331              : 
     332              :             CALL tfw_p_3(rho(ispin)%array, norm_drho(ispin)%array, &
     333              :                          rho_1_3(ispin)%array, s, e_rho_rho_rho, e_rho_rho_ndrho, &
     334            0 :                          e_rho_ndrho_ndrho, npoints)
     335              :          END IF
     336            0 :          IF (order > 3 .OR. order < -3) THEN
     337            0 :             CPABORT("derivatives bigger than 3 not implemented")
     338              :          END IF
     339              :       END DO
     340              : 
     341            0 :       DEALLOCATE (s)
     342            0 :       CALL timestop(handle)
     343            0 :    END SUBROUTINE tfw_lsd_eval
     344              : 
     345              : ! **************************************************************************************************
     346              : !> \brief ...
     347              : !> \param rho ...
     348              : !> \param r13 ...
     349              : !> \param s ...
     350              : !> \param e_0 ...
     351              : !> \param npoints ...
     352              : ! **************************************************************************************************
     353            0 :    SUBROUTINE tfw_u_0(rho, r13, s, e_0, npoints)
     354              : 
     355              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13, s
     356              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0
     357              :       INTEGER, INTENT(in)                                :: npoints
     358              : 
     359              :       INTEGER                                            :: ip
     360              : 
     361              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     362            0 : !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13,s,fvw)
     363              :       DO ip = 1, npoints
     364              : 
     365              :          IF (rho(ip) > eps_rho) THEN
     366              : 
     367              :             e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip) + fvw*s(ip)
     368              : 
     369              :          END IF
     370              : 
     371              :       END DO
     372              : 
     373            0 :    END SUBROUTINE tfw_u_0
     374              : 
     375              : ! **************************************************************************************************
     376              : !> \brief ...
     377              : !> \param rho ...
     378              : !> \param grho ...
     379              : !> \param r13 ...
     380              : !> \param s ...
     381              : !> \param e_rho ...
     382              : !> \param e_ndrho ...
     383              : !> \param npoints ...
     384              : ! **************************************************************************************************
     385            0 :    SUBROUTINE tfw_u_1(rho, grho, r13, s, e_rho, e_ndrho, npoints)
     386              : 
     387              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, grho, r13, s
     388              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho, e_ndrho
     389              :       INTEGER, INTENT(in)                                :: npoints
     390              : 
     391              :       INTEGER                                            :: ip
     392              :       REAL(KIND=dp)                                      :: f
     393              : 
     394            0 :       f = f53*flda
     395              : 
     396              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     397            0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho,e_ndrho,grho,s,r13,f,fvw)
     398              :       DO ip = 1, npoints
     399              : 
     400              :          IF (rho(ip) > eps_rho) THEN
     401              : 
     402              :             e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip) - fvw*s(ip)/rho(ip)
     403              :             e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grho(ip)/rho(ip)
     404              : 
     405              :          END IF
     406              : 
     407              :       END DO
     408              : 
     409            0 :    END SUBROUTINE tfw_u_1
     410              : 
     411              : ! **************************************************************************************************
     412              : !> \brief ...
     413              : !> \param rho ...
     414              : !> \param grho ...
     415              : !> \param r13 ...
     416              : !> \param s ...
     417              : !> \param e_rho_rho ...
     418              : !> \param e_rho_ndrho ...
     419              : !> \param e_ndrho_ndrho ...
     420              : !> \param npoints ...
     421              : ! **************************************************************************************************
     422            0 :    SUBROUTINE tfw_u_2(rho, grho, r13, s, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
     423              :                       npoints)
     424              : 
     425              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, grho, r13, s
     426              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
     427              :       INTEGER, INTENT(in)                                :: npoints
     428              : 
     429              :       INTEGER                                            :: ip
     430              :       REAL(KIND=dp)                                      :: f
     431              : 
     432            0 :       f = f23*f53*flda
     433              : 
     434              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     435            0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho,grho,f,fvw)
     436              :       DO ip = 1, npoints
     437              : 
     438              :          IF (rho(ip) > eps_rho) THEN
     439              : 
     440              :             e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip) + 2.0_dp*fvw*s(ip)/(rho(ip)*rho(ip))
     441              :             e_rho_ndrho(ip) = e_rho_ndrho(ip) - 2.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip))
     442              :             e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rho(ip)
     443              : 
     444              :          END IF
     445              : 
     446              :       END DO
     447              : 
     448            0 :    END SUBROUTINE tfw_u_2
     449              : 
     450              : ! **************************************************************************************************
     451              : !> \brief ...
     452              : !> \param rho ...
     453              : !> \param grho ...
     454              : !> \param r13 ...
     455              : !> \param s ...
     456              : !> \param e_rho_rho_rho ...
     457              : !> \param e_rho_rho_ndrho ...
     458              : !> \param e_rho_ndrho_ndrho ...
     459              : !> \param npoints ...
     460              : ! **************************************************************************************************
     461            0 :    SUBROUTINE tfw_u_3(rho, grho, r13, s, e_rho_rho_rho, e_rho_rho_ndrho, &
     462              :                       e_rho_ndrho_ndrho, npoints)
     463              : 
     464              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, grho, r13, s
     465              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho_rho, e_rho_rho_ndrho, &
     466              :                                                             e_rho_ndrho_ndrho
     467              :       INTEGER, INTENT(in)                                :: npoints
     468              : 
     469              :       INTEGER                                            :: ip
     470              :       REAL(KIND=dp)                                      :: f
     471              : 
     472            0 :       f = -f13*f23*f53*flda
     473              : 
     474              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     475            0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,r13,s,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw)
     476              :       DO ip = 1, npoints
     477              : 
     478              :          IF (rho(ip) > eps_rho) THEN
     479              : 
     480              :             e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip)) &
     481              :                                 - 6.0_dp*fvw*s(ip)/(rho(ip)*rho(ip)*rho(ip))
     482              :             e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
     483              :                                   + 4.0_dp*fvw*grho(ip)/(rho(ip)*rho(ip)*rho(ip))
     484              :             e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
     485              :                                     - 2.0_dp*fvw/(rho(ip)*rho(ip))
     486              :          END IF
     487              : 
     488              :       END DO
     489              : 
     490            0 :    END SUBROUTINE tfw_u_3
     491              : 
     492              : ! **************************************************************************************************
     493              : !> \brief ...
     494              : !> \param rhoa ...
     495              : !> \param r13a ...
     496              : !> \param sa ...
     497              : !> \param e_0 ...
     498              : !> \param npoints ...
     499              : ! **************************************************************************************************
     500            0 :    SUBROUTINE tfw_p_0(rhoa, r13a, sa, e_0, npoints)
     501              : 
     502              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, r13a, sa
     503              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0
     504              :       INTEGER, INTENT(in)                                :: npoints
     505              : 
     506              :       INTEGER                                            :: ip
     507              : 
     508              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     509            0 : !$OMP SHARED(npoints, rhoa,eps_rho,e_0,r13a,sa,flsd,fvw)
     510              :       DO ip = 1, npoints
     511              : 
     512              :          IF (rhoa(ip) > eps_rho) THEN
     513              :             e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip) + fvw*sa(ip)
     514              :          END IF
     515              : 
     516              :       END DO
     517              : 
     518            0 :    END SUBROUTINE tfw_p_0
     519              : 
     520              : ! **************************************************************************************************
     521              : !> \brief ...
     522              : !> \param rhoa ...
     523              : !> \param grhoa ...
     524              : !> \param r13a ...
     525              : !> \param sa ...
     526              : !> \param e_rho ...
     527              : !> \param e_ndrho ...
     528              : !> \param npoints ...
     529              : ! **************************************************************************************************
     530            0 :    SUBROUTINE tfw_p_1(rhoa, grhoa, r13a, sa, e_rho, e_ndrho, npoints)
     531              : 
     532              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, grhoa, r13a, sa
     533              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho, e_ndrho
     534              :       INTEGER, INTENT(in)                                :: npoints
     535              : 
     536              :       INTEGER                                            :: ip
     537              :       REAL(KIND=dp)                                      :: f
     538              : 
     539            0 :       f = f53*flsd
     540              : 
     541              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     542            0 : !$OMP SHARED(npoints,rhoa,eps_rho,r13a,sa,fvw,grhoa,e_rho,e_ndrho,f)
     543              :       DO ip = 1, npoints
     544              : 
     545              :          IF (rhoa(ip) > eps_rho) THEN
     546              :             e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip) - fvw*sa(ip)/rhoa(ip)
     547              :             e_ndrho(ip) = e_ndrho(ip) + 2.0_dp*fvw*grhoa(ip)/rhoa(ip)
     548              :          END IF
     549              : 
     550              :       END DO
     551              : 
     552            0 :    END SUBROUTINE tfw_p_1
     553              : 
     554              : ! **************************************************************************************************
     555              : !> \brief ...
     556              : !> \param rhoa ...
     557              : !> \param grhoa ...
     558              : !> \param r13a ...
     559              : !> \param sa ...
     560              : !> \param e_rho_rho ...
     561              : !> \param e_rho_ndrho ...
     562              : !> \param e_ndrho_ndrho ...
     563              : !> \param npoints ...
     564              : ! **************************************************************************************************
     565            0 :    SUBROUTINE tfw_p_2(rhoa, grhoa, r13a, sa, e_rho_rho, e_rho_ndrho, &
     566              :                       e_ndrho_ndrho, npoints)
     567              : 
     568              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, grhoa, r13a, sa
     569              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
     570              :       INTEGER, INTENT(in)                                :: npoints
     571              : 
     572              :       INTEGER                                            :: ip
     573              :       REAL(KIND=dp)                                      :: f
     574              : 
     575            0 :       f = f23*f53*flsd
     576              : 
     577              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     578            0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,fvw,r13a,sa,e_rho_ndrho,e_ndrho_ndrho)
     579              :       DO ip = 1, npoints
     580              : 
     581              :          IF (rhoa(ip) > eps_rho) THEN
     582              :             e_rho_rho(ip) = e_rho_rho(ip) &
     583              :                             + f/r13a(ip) + 2.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip))
     584              :             e_rho_ndrho(ip) = e_rho_ndrho(ip) &
     585              :                               - 2.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip))
     586              :             e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + 2.0_dp*fvw/rhoa(ip)
     587              :          END IF
     588              : 
     589              :       END DO
     590              : 
     591            0 :    END SUBROUTINE tfw_p_2
     592              : 
     593              : ! **************************************************************************************************
     594              : !> \brief ...
     595              : !> \param rhoa ...
     596              : !> \param grhoa ...
     597              : !> \param r13a ...
     598              : !> \param sa ...
     599              : !> \param e_rho_rho_rho ...
     600              : !> \param e_rho_rho_ndrho ...
     601              : !> \param e_rho_ndrho_ndrho ...
     602              : !> \param npoints ...
     603              : ! **************************************************************************************************
     604            0 :    SUBROUTINE tfw_p_3(rhoa, grhoa, r13a, sa, e_rho_rho_rho, e_rho_rho_ndrho, &
     605              :                       e_rho_ndrho_ndrho, npoints)
     606              : 
     607              :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, grhoa, r13a, sa
     608              :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho_rho, e_rho_rho_ndrho, &
     609              :                                                             e_rho_ndrho_ndrho
     610              :       INTEGER, INTENT(in)                                :: npoints
     611              : 
     612              :       INTEGER                                            :: ip
     613              :       REAL(KIND=dp)                                      :: f
     614              : 
     615            0 :       f = -f13*f23*f53*flsd
     616              : 
     617              : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     618            0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,f,fvw,sa,grhoa)
     619              :       DO ip = 1, npoints
     620              : 
     621              :          IF (rhoa(ip) > eps_rho) THEN
     622              :             e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
     623              :                                 + f/(r13a(ip)*rhoa(ip)) &
     624              :                                 - 6.0_dp*fvw*sa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
     625              :             e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
     626              :                                   + 4.0_dp*fvw*grhoa(ip)/(rhoa(ip)*rhoa(ip)*rhoa(ip))
     627              :             e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
     628              :                                     - 2.0_dp*fvw/(rhoa(ip)*rhoa(ip))
     629              :          END IF
     630              : 
     631              :       END DO
     632              : 
     633            0 :    END SUBROUTINE tfw_p_3
     634              : 
     635              : END MODULE xc_tfw
     636              : 
        

Generated by: LCOV version 2.0-1