LCOV - code coverage report
Current view: top level - src/xc - xc_thomas_fermi.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 43 114 37.7 %
Date: 2024-03-29 07:50:05 Functions: 5 12 41.7 %

          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 Calculate the Thomas-Fermi kinetic energy functional
      10             : !> \note
      11             : !>      Order of derivatives is: LDA 0; 1; 2; 3;
      12             : !>                               LSD 0; a  b; aa bb; aaa bbb;
      13             : !> \par History
      14             : !>      JGH (26.02.2003) : OpenMP enabled
      15             : !>      fawzi (04.2004)  : adapted to the new xc interface
      16             : !> \author JGH (18.02.2002)
      17             : ! **************************************************************************************************
      18             : MODULE xc_thomas_fermi
      19             :    USE cp_array_utils,                  ONLY: cp_3d_r_cp_type
      20             :    USE kinds,                           ONLY: dp
      21             :    USE xc_derivative_desc,              ONLY: deriv_rho,&
      22             :                                               deriv_rhoa,&
      23             :                                               deriv_rhob
      24             :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      25             :                                               xc_dset_get_derivative
      26             :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      27             :                                               xc_derivative_type
      28             :    USE xc_functionals_utilities,        ONLY: set_util
      29             :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      30             :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      31             :                                               xc_rho_set_type
      32             : #include "../base/base_uses.f90"
      33             : 
      34             :    IMPLICIT NONE
      35             : 
      36             :    PRIVATE
      37             : 
      38             : ! *** Global parameters ***
      39             : 
      40             :    REAL(KIND=dp), PARAMETER :: pi = 3.14159265358979323846264338_dp
      41             :    REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
      42             :                                f23 = 2.0_dp*f13, &
      43             :                                f43 = 4.0_dp*f13, &
      44             :                                f53 = 5.0_dp*f13
      45             : 
      46             :    PUBLIC :: thomas_fermi_info, thomas_fermi_lda_eval, thomas_fermi_lsd_eval
      47             : 
      48             :    REAL(KIND=dp) :: cf, flda, flsd
      49             :    REAL(KIND=dp) :: eps_rho
      50             : 
      51             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_thomas_fermi'
      52             : 
      53             : CONTAINS
      54             : 
      55             : ! **************************************************************************************************
      56             : !> \brief ...
      57             : !> \param cutoff ...
      58             : ! **************************************************************************************************
      59         216 :    SUBROUTINE thomas_fermi_init(cutoff)
      60             : 
      61             :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
      62             : 
      63         216 :       eps_rho = cutoff
      64         216 :       CALL set_util(cutoff)
      65             : 
      66         216 :       cf = 0.3_dp*(3.0_dp*pi*pi)**f23
      67         216 :       flda = cf
      68         216 :       flsd = flda*2.0_dp**f23
      69             : 
      70         216 :    END SUBROUTINE thomas_fermi_init
      71             : 
      72             : ! **************************************************************************************************
      73             : !> \brief ...
      74             : !> \param lsd ...
      75             : !> \param reference ...
      76             : !> \param shortform ...
      77             : !> \param needs ...
      78             : !> \param max_deriv ...
      79             : ! **************************************************************************************************
      80         216 :    SUBROUTINE thomas_fermi_info(lsd, reference, shortform, needs, max_deriv)
      81             :       LOGICAL, INTENT(in)                                :: lsd
      82             :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      83             :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      84             :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      85             : 
      86         216 :       IF (PRESENT(reference)) THEN
      87           0 :          reference = "Thomas-Fermi kinetic energy functional: see Parr and Yang"
      88           0 :          IF (.NOT. lsd) THEN
      89           0 :             IF (LEN_TRIM(reference) + 6 < LEN(reference)) THEN
      90           0 :                reference(LEN_TRIM(reference):LEN_TRIM(reference) + 6) = ' {LDA}'
      91             :             END IF
      92             :          END IF
      93             :       END IF
      94         216 :       IF (PRESENT(shortform)) THEN
      95           0 :          shortform = "Thomas-Fermi kinetic energy functional"
      96           0 :          IF (.NOT. lsd) THEN
      97           0 :             IF (LEN_TRIM(shortform) + 6 < LEN(shortform)) THEN
      98           0 :                shortform(LEN_TRIM(shortform):LEN_TRIM(shortform) + 6) = ' {LDA}'
      99             :             END IF
     100             :          END IF
     101             :       END IF
     102         216 :       IF (PRESENT(needs)) THEN
     103         216 :          IF (lsd) THEN
     104           0 :             needs%rho_spin = .TRUE.
     105           0 :             needs%rho_spin_1_3 = .TRUE.
     106             :          ELSE
     107         216 :             needs%rho = .TRUE.
     108         216 :             needs%rho_1_3 = .TRUE.
     109             :          END IF
     110             :       END IF
     111         216 :       IF (PRESENT(max_deriv)) max_deriv = 3
     112             : 
     113         216 :    END SUBROUTINE thomas_fermi_info
     114             : 
     115             : ! **************************************************************************************************
     116             : !> \brief ...
     117             : !> \param rho_set ...
     118             : !> \param deriv_set ...
     119             : !> \param order ...
     120             : ! **************************************************************************************************
     121         432 :    SUBROUTINE thomas_fermi_lda_eval(rho_set, deriv_set, order)
     122             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     123             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     124             :       INTEGER, INTENT(in)                                :: order
     125             : 
     126             :       CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lda_eval'
     127             : 
     128             :       INTEGER                                            :: handle, npoints
     129             :       INTEGER, DIMENSION(2, 3)                           :: bo
     130             :       REAL(KIND=dp)                                      :: epsilon_rho
     131             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     132         216 :          POINTER                                         :: e_0, e_rho, e_rho_rho, e_rho_rho_rho, &
     133         216 :                                                             r13, rho
     134             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     135             : 
     136         216 :       CALL timeset(routineN, handle)
     137             : 
     138             :       CALL xc_rho_set_get(rho_set, rho_1_3=r13, rho=rho, &
     139         216 :                           local_bounds=bo, rho_cutoff=epsilon_rho)
     140         216 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     141         216 :       CALL thomas_fermi_init(epsilon_rho)
     142             : 
     143         216 :       IF (order >= 0) THEN
     144             :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     145         216 :                                          allocate_deriv=.TRUE.)
     146         216 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     147             : 
     148         216 :          CALL thomas_fermi_lda_0(rho, r13, e_0, npoints)
     149             :       END IF
     150         216 :       IF (order >= 1 .OR. order == -1) THEN
     151             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     152         216 :                                          allocate_deriv=.TRUE.)
     153         216 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     154             : 
     155         216 :          CALL thomas_fermi_lda_1(rho, r13, e_rho, npoints)
     156             :       END IF
     157         216 :       IF (order >= 2 .OR. order == -2) THEN
     158             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
     159           0 :                                          allocate_deriv=.TRUE.)
     160           0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     161             : 
     162           0 :          CALL thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
     163             :       END IF
     164         216 :       IF (order >= 3 .OR. order == -3) THEN
     165             :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
     166           0 :                                          allocate_deriv=.TRUE.)
     167           0 :          CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     168             : 
     169           0 :          CALL thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
     170             :       END IF
     171         216 :       IF (order > 3 .OR. order < -3) THEN
     172           0 :          CPABORT("derivatives bigger than 3 not implemented")
     173             :       END IF
     174         216 :       CALL timestop(handle)
     175         216 :    END SUBROUTINE thomas_fermi_lda_eval
     176             : 
     177             : ! **************************************************************************************************
     178             : !> \brief ...
     179             : !> \param rho_set ...
     180             : !> \param deriv_set ...
     181             : !> \param order ...
     182             : ! **************************************************************************************************
     183           0 :    SUBROUTINE thomas_fermi_lsd_eval(rho_set, deriv_set, order)
     184             :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     185             :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     186             :       INTEGER, INTENT(in)                                :: order
     187             : 
     188             :       CHARACTER(len=*), PARAMETER :: routineN = 'thomas_fermi_lsd_eval'
     189             :       INTEGER, DIMENSION(2), PARAMETER :: rho_spin_name = [deriv_rhoa, deriv_rhob]
     190             : 
     191             :       INTEGER                                            :: handle, i, ispin, npoints
     192             :       INTEGER, DIMENSION(2, 3)                           :: bo
     193             :       REAL(KIND=dp)                                      :: epsilon_rho
     194             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     195           0 :          POINTER                                         :: e_0, e_rho, e_rho_rho, e_rho_rho_rho
     196           0 :       TYPE(cp_3d_r_cp_type), DIMENSION(2)                :: rho, rho_1_3
     197             :       TYPE(xc_derivative_type), POINTER                  :: deriv
     198             : 
     199           0 :       CALL timeset(routineN, handle)
     200           0 :       NULLIFY (deriv)
     201           0 :       DO i = 1, 2
     202           0 :          NULLIFY (rho(i)%array, rho_1_3(i)%array)
     203             :       END DO
     204             : 
     205             :       CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
     206             :                           rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
     207             :                           rhob=rho(2)%array, &
     208             :                           rho_cutoff=epsilon_rho, &
     209           0 :                           local_bounds=bo)
     210           0 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     211           0 :       CALL thomas_fermi_init(epsilon_rho)
     212             : 
     213           0 :       DO ispin = 1, 2
     214           0 :          IF (order >= 0) THEN
     215             :             deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     216           0 :                                             allocate_deriv=.TRUE.)
     217           0 :             CALL xc_derivative_get(deriv, deriv_data=e_0)
     218             : 
     219             :             CALL thomas_fermi_lsd_0(rho(ispin)%array, rho_1_3(ispin)%array, &
     220           0 :                                     e_0, npoints)
     221             :          END IF
     222           0 :          IF (order >= 1 .OR. order == -1) THEN
     223             :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin)], &
     224           0 :                                             allocate_deriv=.TRUE.)
     225           0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho)
     226             : 
     227             :             CALL thomas_fermi_lsd_1(rho(ispin)%array, rho_1_3(ispin)%array, &
     228           0 :                                     e_rho, npoints)
     229             :          END IF
     230           0 :          IF (order >= 2 .OR. order == -2) THEN
     231             :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     232           0 :                                                         rho_spin_name(ispin)], allocate_deriv=.TRUE.)
     233           0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
     234             : 
     235             :             CALL thomas_fermi_lsd_2(rho(ispin)%array, rho_1_3(ispin)%array, &
     236           0 :                                     e_rho_rho, npoints)
     237             :          END IF
     238           0 :          IF (order >= 3 .OR. order == -3) THEN
     239             :             deriv => xc_dset_get_derivative(deriv_set, [rho_spin_name(ispin), &
     240             :                                                         rho_spin_name(ispin), rho_spin_name(ispin)], &
     241           0 :                                             allocate_deriv=.TRUE.)
     242           0 :             CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
     243             : 
     244             :             CALL thomas_fermi_lsd_3(rho(ispin)%array, rho_1_3(ispin)%array, &
     245           0 :                                     e_rho_rho_rho, npoints)
     246             :          END IF
     247           0 :          IF (order > 3 .OR. order < -3) THEN
     248           0 :             CPABORT("derivatives bigger than 3 not implemented")
     249             :          END IF
     250             :       END DO
     251           0 :       CALL timestop(handle)
     252           0 :    END SUBROUTINE thomas_fermi_lsd_eval
     253             : 
     254             : ! **************************************************************************************************
     255             : !> \brief ...
     256             : !> \param rho ...
     257             : !> \param r13 ...
     258             : !> \param e_0 ...
     259             : !> \param npoints ...
     260             : ! **************************************************************************************************
     261         216 :    SUBROUTINE thomas_fermi_lda_0(rho, r13, e_0, npoints)
     262             : 
     263             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13
     264             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0
     265             :       INTEGER, INTENT(in)                                :: npoints
     266             : 
     267             :       INTEGER                                            :: ip
     268             : 
     269             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     270         216 : !$OMP SHARED(npoints,rho,eps_rho,e_0,flda,r13)
     271             :       DO ip = 1, npoints
     272             : 
     273             :          IF (rho(ip) > eps_rho) THEN
     274             : 
     275             :             e_0(ip) = e_0(ip) + flda*r13(ip)*r13(ip)*rho(ip)
     276             : 
     277             :          END IF
     278             : 
     279             :       END DO
     280             : 
     281         216 :    END SUBROUTINE thomas_fermi_lda_0
     282             : 
     283             : ! **************************************************************************************************
     284             : !> \brief ...
     285             : !> \param rho ...
     286             : !> \param r13 ...
     287             : !> \param e_rho ...
     288             : !> \param npoints ...
     289             : ! **************************************************************************************************
     290         216 :    SUBROUTINE thomas_fermi_lda_1(rho, r13, e_rho, npoints)
     291             : 
     292             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13
     293             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho
     294             :       INTEGER, INTENT(in)                                :: npoints
     295             : 
     296             :       INTEGER                                            :: ip
     297             :       REAL(KIND=dp)                                      :: f
     298             : 
     299         216 :       f = f53*flda
     300             : 
     301             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE) &
     302         216 : !$OMP SHARED(npoints,rho,eps_rho,e_rho,f,r13)
     303             :       DO ip = 1, npoints
     304             : 
     305             :          IF (rho(ip) > eps_rho) THEN
     306             : 
     307             :             e_rho(ip) = e_rho(ip) + f*r13(ip)*r13(ip)
     308             : 
     309             :          END IF
     310             : 
     311             :       END DO
     312             : 
     313         216 :    END SUBROUTINE thomas_fermi_lda_1
     314             : 
     315             : ! **************************************************************************************************
     316             : !> \brief ...
     317             : !> \param rho ...
     318             : !> \param r13 ...
     319             : !> \param e_rho_rho ...
     320             : !> \param npoints ...
     321             : ! **************************************************************************************************
     322           0 :    SUBROUTINE thomas_fermi_lda_2(rho, r13, e_rho_rho, npoints)
     323             : 
     324             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13
     325             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho
     326             :       INTEGER, INTENT(in)                                :: npoints
     327             : 
     328             :       INTEGER                                            :: ip
     329             :       REAL(KIND=dp)                                      :: f
     330             : 
     331           0 :       f = f23*f53*flda
     332             : 
     333             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     334           0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho,f,r13)
     335             :       DO ip = 1, npoints
     336             : 
     337             :          IF (rho(ip) > eps_rho) THEN
     338             : 
     339             :             e_rho_rho(ip) = e_rho_rho(ip) + f/r13(ip)
     340             : 
     341             :          END IF
     342             : 
     343             :       END DO
     344             : 
     345           0 :    END SUBROUTINE thomas_fermi_lda_2
     346             : 
     347             : ! **************************************************************************************************
     348             : !> \brief ...
     349             : !> \param rho ...
     350             : !> \param r13 ...
     351             : !> \param e_rho_rho_rho ...
     352             : !> \param npoints ...
     353             : ! **************************************************************************************************
     354           0 :    SUBROUTINE thomas_fermi_lda_3(rho, r13, e_rho_rho_rho, npoints)
     355             : 
     356             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rho, r13
     357             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho_rho
     358             :       INTEGER, INTENT(in)                                :: npoints
     359             : 
     360             :       INTEGER                                            :: ip
     361             :       REAL(KIND=dp)                                      :: f
     362             : 
     363           0 :       f = -f13*f23*f53*flda
     364             : 
     365             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     366           0 : !$OMP SHARED(npoints,rho,eps_rho,e_rho_rho_rho,f,r13)
     367             :       DO ip = 1, npoints
     368             : 
     369             :          IF (rho(ip) > eps_rho) THEN
     370             : 
     371             :             e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13(ip)*rho(ip))
     372             : 
     373             :          END IF
     374             : 
     375             :       END DO
     376             : 
     377           0 :    END SUBROUTINE thomas_fermi_lda_3
     378             : 
     379             : ! **************************************************************************************************
     380             : !> \brief ...
     381             : !> \param rhoa ...
     382             : !> \param r13a ...
     383             : !> \param e_0 ...
     384             : !> \param npoints ...
     385             : ! **************************************************************************************************
     386           0 :    SUBROUTINE thomas_fermi_lsd_0(rhoa, r13a, e_0, npoints)
     387             : 
     388             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, r13a
     389             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_0
     390             :       INTEGER, INTENT(in)                                :: npoints
     391             : 
     392             :       INTEGER                                            :: ip
     393             : 
     394             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     395           0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_0,flsd,r13a)
     396             :       DO ip = 1, npoints
     397             : 
     398             :          IF (rhoa(ip) > eps_rho) THEN
     399             :             e_0(ip) = e_0(ip) + flsd*r13a(ip)*r13a(ip)*rhoa(ip)
     400             :          END IF
     401             : 
     402             :       END DO
     403             : 
     404           0 :    END SUBROUTINE thomas_fermi_lsd_0
     405             : 
     406             : ! **************************************************************************************************
     407             : !> \brief ...
     408             : !> \param rhoa ...
     409             : !> \param r13a ...
     410             : !> \param e_rho ...
     411             : !> \param npoints ...
     412             : ! **************************************************************************************************
     413           0 :    SUBROUTINE thomas_fermi_lsd_1(rhoa, r13a, e_rho, npoints)
     414             : 
     415             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, r13a
     416             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho
     417             :       INTEGER, INTENT(in)                                :: npoints
     418             : 
     419             :       INTEGER                                            :: ip
     420             :       REAL(KIND=dp)                                      :: f
     421             : 
     422           0 :       f = f53*flsd
     423             : 
     424             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     425           0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho,f,r13a)
     426             :       DO ip = 1, npoints
     427             : 
     428             :          IF (rhoa(ip) > eps_rho) THEN
     429             :             e_rho(ip) = e_rho(ip) + f*r13a(ip)*r13a(ip)
     430             :          END IF
     431             : 
     432             :       END DO
     433             : 
     434           0 :    END SUBROUTINE thomas_fermi_lsd_1
     435             : 
     436             : ! **************************************************************************************************
     437             : !> \brief ...
     438             : !> \param rhoa ...
     439             : !> \param r13a ...
     440             : !> \param e_rho_rho ...
     441             : !> \param npoints ...
     442             : ! **************************************************************************************************
     443           0 :    SUBROUTINE thomas_fermi_lsd_2(rhoa, r13a, e_rho_rho, npoints)
     444             : 
     445             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, r13a
     446             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho
     447             :       INTEGER, INTENT(in)                                :: npoints
     448             : 
     449             :       INTEGER                                            :: ip
     450             :       REAL(KIND=dp)                                      :: f
     451             : 
     452           0 :       f = f23*f53*flsd
     453             : 
     454             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     455           0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho,f,r13a)
     456             : 
     457             :       DO ip = 1, npoints
     458             : 
     459             :          IF (rhoa(ip) > eps_rho) THEN
     460             :             e_rho_rho(ip) = e_rho_rho(ip) + f/r13a(ip)
     461             :          END IF
     462             : 
     463             :       END DO
     464             : 
     465           0 :    END SUBROUTINE thomas_fermi_lsd_2
     466             : 
     467             : ! **************************************************************************************************
     468             : !> \brief ...
     469             : !> \param rhoa ...
     470             : !> \param r13a ...
     471             : !> \param e_rho_rho_rho ...
     472             : !> \param npoints ...
     473             : ! **************************************************************************************************
     474           0 :    SUBROUTINE thomas_fermi_lsd_3(rhoa, r13a, e_rho_rho_rho, npoints)
     475             : 
     476             :       REAL(KIND=dp), DIMENSION(*), INTENT(IN)            :: rhoa, r13a
     477             :       REAL(KIND=dp), DIMENSION(*), INTENT(INOUT)         :: e_rho_rho_rho
     478             :       INTEGER, INTENT(in)                                :: npoints
     479             : 
     480             :       INTEGER                                            :: ip
     481             :       REAL(KIND=dp)                                      :: f
     482             : 
     483           0 :       f = -f13*f23*f53*flsd
     484             : 
     485             : !$OMP PARALLEL DO PRIVATE(ip) DEFAULT(NONE)&
     486           0 : !$OMP SHARED(npoints,rhoa,eps_rho,e_rho_rho_rho,f,r13a)
     487             :       DO ip = 1, npoints
     488             : 
     489             :          IF (rhoa(ip) > eps_rho) THEN
     490             :             e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + f/(r13a(ip)*rhoa(ip))
     491             :          END IF
     492             : 
     493             :       END DO
     494             : 
     495           0 :    END SUBROUTINE thomas_fermi_lsd_3
     496             : 
     497             : END MODULE xc_thomas_fermi
     498             : 

Generated by: LCOV version 1.15