LCOV - code coverage report
Current view: top level - src/xc - xc_tpss.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.7 % 377 372
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            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 Calculates the tpss functional.
      10              : !> \note
      11              : !>      The derivation of the formulaes is lengthly, and not fully trivial,
      12              : !>      so I have put it in doc/tpss.mw
      13              : !> \par History
      14              : !>      05.2004 created
      15              : !> \author fawzi
      16              : ! **************************************************************************************************
      17              : MODULE xc_tpss
      18              :    USE bibliography,                    ONLY: Tao2003,&
      19              :                                               cite_reference
      20              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      21              :                                               cp_logger_type
      22              :    USE input_section_types,             ONLY: section_vals_type,&
      23              :                                               section_vals_val_get
      24              :    USE kinds,                           ONLY: dp
      25              :    USE mathconstants,                   ONLY: pi
      26              :    USE xc_derivative_desc,              ONLY: deriv_norm_drho,&
      27              :                                               deriv_rho,&
      28              :                                               deriv_tau
      29              :    USE xc_derivative_set_types,         ONLY: xc_derivative_set_type,&
      30              :                                               xc_dset_get_derivative
      31              :    USE xc_derivative_types,             ONLY: xc_derivative_get,&
      32              :                                               xc_derivative_type
      33              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      34              :    USE xc_rho_set_types,                ONLY: xc_rho_set_get,&
      35              :                                               xc_rho_set_type
      36              : #include "../base/base_uses.f90"
      37              : 
      38              :    IMPLICIT NONE
      39              :    PRIVATE
      40              : 
      41              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tpss'
      43              : 
      44              :    PUBLIC :: tpss_lda_info, tpss_lda_eval
      45              : 
      46              : !***
      47              : CONTAINS
      48              : 
      49              : ! **************************************************************************************************
      50              : !> \brief return various information on the functional
      51              : !> \param tpss_params ...
      52              : !> \param reference string with the reference of the actual functional
      53              : !> \param shortform string with the shortform of the functional name
      54              : !> \param needs the components needed by this functional are set to
      55              : !>        true (does not set the unneeded components to false)
      56              : !> \param max_deriv the highest derivative available
      57              : !> \author fawzi
      58              : ! **************************************************************************************************
      59         1692 :    SUBROUTINE tpss_lda_info(tpss_params, reference, shortform, needs, max_deriv)
      60              :       TYPE(section_vals_type), POINTER                   :: tpss_params
      61              :       CHARACTER(LEN=*), INTENT(OUT), OPTIONAL            :: reference, shortform
      62              :       TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL  :: needs
      63              :       INTEGER, INTENT(out), OPTIONAL                     :: max_deriv
      64              : 
      65              :       REAL(kind=dp)                                      :: sc, sx
      66              : 
      67          846 :       CALL section_vals_val_get(tpss_params, "SCALE_C", r_val=sc)
      68          846 :       CALL section_vals_val_get(tpss_params, "SCALE_X", r_val=sx)
      69              : 
      70          846 :       IF (PRESENT(reference)) THEN
      71           12 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
      72           12 :             reference = "J. Tao, J.P.Perdew, V.N.Staroverov, E.Scuseria PRL, 91, 146401 (2003) {LDA version}"
      73              :          ELSE
      74              :             WRITE (reference, "(a,'sx=',f5.3,'sc=',f5.3,' {LDA version}')") &
      75            0 :                "J. Tao, J.P.Perdew, V.N.Staroverov, E.Scuseria PRL, 91, 146401 (2003)", &
      76            0 :                sx, sc
      77              :          END IF
      78              :       END IF
      79          846 :       IF (PRESENT(shortform)) THEN
      80           12 :          IF (sx == 1._dp .AND. sc == 1._dp) THEN
      81           12 :             shortform = "TPSS meta-GGA functional (LDA)"
      82              :          ELSE
      83              :             WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,' (LDA)')") &
      84            0 :                "TPSS meta-GGA functional", &
      85            0 :                sx, sc
      86              :          END IF
      87              :       END IF
      88          846 :       IF (PRESENT(needs)) THEN
      89          834 :          needs%rho = .TRUE.
      90          834 :          needs%tau = .TRUE.
      91          834 :          needs%norm_drho = .TRUE.
      92              :       END IF
      93          846 :       IF (PRESENT(max_deriv)) max_deriv = 1
      94              : 
      95          846 :    END SUBROUTINE tpss_lda_info
      96              : 
      97              : ! **************************************************************************************************
      98              : !> \brief evaluates the tpss functional in the spin unpolarized (lda) case
      99              : !> \param rho_set the density where you want to evaluate the functional
     100              : !> \param deriv_set place where to store the functional derivatives (they are
     101              : !>        added to the derivatives)
     102              : !> \param grad_deriv degree of the derivative that should be evaluated,
     103              : !>        if positive all the derivatives up to the given degree are evaluated,
     104              : !>        if negative only the given degree is calculated
     105              : !> \param tpss_params ...
     106              : !> \author fawzi
     107              : ! **************************************************************************************************
     108         4248 :    SUBROUTINE tpss_lda_eval(rho_set, deriv_set, grad_deriv, tpss_params)
     109              :       TYPE(xc_rho_set_type), INTENT(IN)                  :: rho_set
     110              :       TYPE(xc_derivative_set_type), INTENT(IN)           :: deriv_set
     111              :       INTEGER, INTENT(in)                                :: grad_deriv
     112              :       TYPE(section_vals_type), POINTER                   :: tpss_params
     113              : 
     114              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'tpss_lda_eval'
     115              : 
     116              :       INTEGER                                            :: handle, non_coer, npoints
     117              :       INTEGER, DIMENSION(2, 3)                           :: bo
     118              :       REAL(kind=dp)                                      :: epsilon_rho, epsilon_tau, scale_ec, &
     119              :                                                             scale_ex
     120              :       REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
     121         1062 :          POINTER                                         :: dummy, e_0, e_ndrho, e_rho, e_tau, &
     122         1062 :                                                             norm_drho, rho, tau
     123              :       TYPE(cp_logger_type), POINTER                      :: logger
     124              :       TYPE(xc_derivative_type), POINTER                  :: deriv
     125              : 
     126         1062 :       CALL timeset(routineN, handle)
     127              : 
     128         1062 :       CALL cite_reference(tao2003)
     129              : 
     130              :       CALL xc_rho_set_get(rho_set, rho=rho, &
     131              :                           norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
     132         1062 :                           tau=tau, tau_cutoff=epsilon_tau)
     133         1062 :       npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
     134              : 
     135         1062 :       dummy => rho
     136              : 
     137         1062 :       e_0 => dummy
     138         1062 :       e_rho => dummy
     139         1062 :       e_ndrho => dummy
     140         1062 :       e_tau => dummy
     141              : 
     142         1062 :       IF (grad_deriv >= 0) THEN
     143              :          deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
     144         1062 :                                          allocate_deriv=.TRUE.)
     145         1062 :          CALL xc_derivative_get(deriv, deriv_data=e_0)
     146              :       END IF
     147         1062 :       IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     148              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
     149         1062 :                                          allocate_deriv=.TRUE.)
     150         1062 :          CALL xc_derivative_get(deriv, deriv_data=e_rho)
     151              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
     152         1062 :                                          allocate_deriv=.TRUE.)
     153         1062 :          CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
     154              :          deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
     155         1062 :                                          allocate_deriv=.TRUE.)
     156         1062 :          CALL xc_derivative_get(deriv, deriv_data=e_tau)
     157              :       END IF
     158         1062 :       IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
     159            0 :          CPABORT("derivatives bigger than 1 not implemented")
     160              :       END IF
     161              : 
     162         1062 :       non_coer = 0
     163         1062 :       CALL section_vals_val_get(tpss_params, "SCALE_C", r_val=scale_ec)
     164         1062 :       CALL section_vals_val_get(tpss_params, "SCALE_X", r_val=scale_ex)
     165              : 
     166              : !$OMP     PARALLEL DEFAULT(NONE) &
     167              : !$OMP              SHARED(rho, tau, norm_drho, e_0, e_rho, e_ndrho, e_tau) &
     168              : !$OMP              SHARED(epsilon_rho, epsilon_tau, npoints, grad_deriv) &
     169              : !$OMP              SHARED(scale_ec, scale_ex) &
     170         1062 : !$OMP              REDUCTION(+: non_coer)
     171              : 
     172              :       CALL tpss_lda_calc(rho=rho, norm_drho=norm_drho, &
     173              :                          tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_tau=e_tau, &
     174              :                          grad_deriv=grad_deriv, npoints=npoints, epsilon_rho=epsilon_rho, &
     175              :                          epsilon_tau=epsilon_tau, scale_ec=scale_ec, scale_ex=scale_ex, non_coer=non_coer)
     176              : 
     177              : !$OMP     END PARALLEL
     178              : 
     179         1062 :       logger => cp_get_default_logger()
     180              :       ! we could check if tau/grad were consistent, but don't do anything here
     181              :       IF (non_coer > 0) THEN
     182         1062 :          non_coer = 0
     183              :       END IF
     184              : 
     185         1062 :       CALL timestop(handle)
     186         1062 :    END SUBROUTINE tpss_lda_eval
     187              : 
     188              : ! **************************************************************************************************
     189              : !> \brief low level calculation routine for the unpolarized (lda) tpss
     190              : !> \param rho ...
     191              : !> \param norm_drho ...
     192              : !> \param tau ...
     193              : !> \param e_0 ...
     194              : !> \param e_rho ...
     195              : !> \param e_ndrho ...
     196              : !> \param e_tau ...
     197              : !> \param npoints ...
     198              : !> \param grad_deriv ...
     199              : !> \param epsilon_rho ...
     200              : !> \param epsilon_tau ...
     201              : !> \param scale_ec ...
     202              : !> \param scale_ex ...
     203              : !> \param non_coer ...
     204              : !> \author fawzi
     205              : !> \note
     206              : !>      maple is nice, but if you want the uman readable version of the code
     207              : !>      look in doc/tpss.mw
     208              : ! **************************************************************************************************
     209         1062 :    SUBROUTINE tpss_lda_calc(rho, norm_drho, tau, e_0, e_rho, e_ndrho, e_tau, &
     210              :                             npoints, grad_deriv, epsilon_rho, epsilon_tau, &
     211              :                             scale_ec, scale_ex, non_coer)
     212              :       REAL(kind=dp), DIMENSION(*), INTENT(in)            :: rho, norm_drho, tau
     213              :       REAL(kind=dp), DIMENSION(*), INTENT(inout)         :: e_0, e_rho, e_ndrho, e_tau
     214              :       INTEGER, INTENT(in)                                :: npoints, grad_deriv
     215              :       REAL(kind=dp), INTENT(in)                          :: epsilon_rho, epsilon_tau, scale_ec, &
     216              :                                                             scale_ex
     217              :       INTEGER, INTENT(inout)                             :: non_coer
     218              : 
     219              :       INTEGER                                            :: abs_grad_deriv, ii
     220              :       LOGICAL                                            :: t571, t639
     221              :       REAL(kind=dp) :: A, A_1, A_2, A_s1, A_s1rho, A_s2, A_s2rho, alpha, alpha_1_1, alpha_1_2, &
     222              :          alphanorm_drho, alpharho, alphatau, Arho, b, beta, beta_1_1, beta_1_2, beta_2_1, &
     223              :          beta_2_2, beta_3_1, beta_3_2, beta_4_1, beta_4_2, c, d, e_c_u_0, e_c_u_0rho, e_c_u_1_s1, &
     224              :          e_c_u_1_s1rho, e_c_u_1_s2, e_c_u_1_s2rho, e_var, epsilon_cGGA, epsilon_cGGA_0_1, &
     225              :          epsilon_cGGA_1_0, epsilon_cGGArho, epsilon_cRevPKZB, epsilon_cRevPKZBnorm_drho, &
     226              :          epsilon_cRevPKZBrho, epsilon_cRevPKZBtau, ex_unif, Fx, gamma_var, Hnorm_drho, k_f_s1, &
     227              :          k_f_s1rho, k_s, k_s_s1, k_s_s2, kappa, m, ma, manorm_drho, marho, mb, mbnorm_drho, mbrho
     228              :       REAL(kind=dp) :: mu, my_ndrho, my_rho, my_tau, p, p_1, p_2, p_3, phi_s1, phi_s2, pnorm_drho, &
     229              :          prho, rs, rs_s1, rs_s1rho, rs_s2, rs_s2rho, rsrho, t, t1, t100, t101, t111, t12, t13, &
     230              :          t138, t14, t140, t143, t145, t146, t147, t151, t152, t16, t161, t168, t177, t186, t187, &
     231              :          t189, t19, t190, t191, t193, t194, t196, t197, t198, t199, t2, t20, t201, t202, t204, &
     232              :          t205, t208, t209, t21, t211, t212, t213, t215, t216, t218, t219, t22, t220, t221, t223, &
     233              :          t224, t226, t227, t230, t231, t233, t234, t235, t238, t239, t241, t242, t243, t245, t246, &
     234              :          t248, t249, t252, t253, t254, t256, t26, t260, t263, t264, t265
     235              :       REAL(kind=dp) :: t267, t268, t269, t27, t271, t272, t274, t275, t276, t277, t278, t279, t28, &
     236              :          t280, t281, t284, t286, t288, t29, t290, t291, t293, t294, t295, t299, t3, t301, t302, &
     237              :          t303, t305, t307, t310, t313, t316, t319, t322, t325, t327, t328, t329, t331, t337, t340, &
     238              :          t343, t344, t35, t351, t36, t370, t371, t376, t383, t385, t386, t39, t390, t391, t395, &
     239              :          t396, t398, t4, t403, t404, t406, t41, t410, t411, t419, t42, t430, t437, t445, t450, &
     240              :          t452, t464, t472, t475, t485, t489, t49, t490, t5, t505, t513, t517, t536, t541, t542, &
     241              :          t546, t547, t549, t55, t554, t555, t557, t561, t562, t569, t574
     242              :       REAL(kind=dp) :: t58, t585, t6, t60, t604, t609, t610, t614, t615, t617, t622, t623, t625, &
     243              :          t629, t630, t637, t642, t645, t659, t67, t7, t71, t73, t77, t78, t79, t799, t80, t84, &
     244              :          t85, t89, t9, t94, t95, t96, t_s1, t_s1norm_drho, t_s1rho, t_s2, t_s2norm_drho, t_s2rho, &
     245              :          tau_w, tau_wnorm_drho, tau_wrho, tildeq_b, tildeq_bnorm_drho, tildeq_brho, tildeq_btau, &
     246              :          tnorm_drho, trho, z, znorm_drho, zrho, ztau
     247              : 
     248              :       IF (.FALSE.) THEN
     249              :          ! useful for testing, we just hack in a well defined functional of tau, ndrho and rho
     250              :          ! and see that things converge properly with OT.
     251              : !$OMP        DO
     252              :          DO ii = 1, npoints
     253              :             my_tau = tau(ii)
     254              :             my_rho = rho(ii)
     255              :             my_ndrho = norm_drho(ii)
     256              :             IF (grad_deriv >= 0) THEN
     257              :                e_0(ii) = e_0(ii) + my_tau*my_ndrho*my_rho
     258              :             END IF
     259              :             IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     260              :                e_rho(ii) = e_rho(ii) + my_tau*my_ndrho
     261              :                e_ndrho(ii) = e_ndrho(ii) + my_tau*my_rho
     262              :                e_tau(ii) = e_tau(ii) + my_rho*my_ndrho
     263              :             END IF
     264              :          END DO
     265              : !$OMP        END DO
     266              :          RETURN
     267              :       END IF
     268              : 
     269         1062 :       abs_grad_deriv = ABS(grad_deriv)
     270              : 
     271         1062 :       kappa = 0.804e0_dp
     272         1062 :       beta = 0.66725e-1_dp
     273         1062 :       mu = 0.21951e0_dp
     274         1062 :       gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2
     275         1062 :       b = 0.4e0_dp
     276         1062 :       c = 0.159096e1_dp
     277         1062 :       e_var = 0.1537e1_dp
     278         1062 :       d = 0.28e1_dp
     279         1062 :       p_1 = 0.10e1_dp
     280         1062 :       A_1 = 0.31091e-1_dp
     281         1062 :       alpha_1_1 = 0.21370e0_dp
     282         1062 :       beta_1_1 = 0.75957e1_dp
     283         1062 :       beta_2_1 = 0.35876e1_dp
     284         1062 :       beta_3_1 = 0.16382e1_dp
     285         1062 :       beta_4_1 = 0.49294e0_dp
     286         1062 :       p_2 = 0.10e1_dp
     287         1062 :       A_2 = 0.15545e-1_dp
     288         1062 :       alpha_1_2 = 0.20548e0_dp
     289         1062 :       beta_1_2 = 0.141189e2_dp
     290         1062 :       beta_2_2 = 0.61977e1_dp
     291         1062 :       beta_3_2 = 0.33662e1_dp
     292         1062 :       beta_4_2 = 0.62517e0_dp
     293         1062 :       p_3 = 0.10e1_dp
     294              : 
     295         1062 :       t1 = 3._dp**(0.1e1_dp/0.3e1_dp)
     296         1062 :       t2 = 4._dp**(0.1e1_dp/0.3e1_dp)
     297         1062 :       t3 = t2**2
     298         1062 :       t4 = t1*t3
     299         1062 :       t5 = 2._dp**(0.1e1_dp/0.3e1_dp)
     300         1062 :       t6 = 0.1e1_dp/pi
     301         1062 :       t12 = t5**2
     302              : 
     303         1062 : !$OMP     DO
     304              : 
     305              :       DO ii = 1, npoints
     306     24682723 :          my_tau = tau(ii)
     307     24682723 :          my_rho = rho(ii)
     308     24682723 :          IF (my_rho > epsilon_rho .AND. my_tau > epsilon_tau) THEN
     309     20673528 :             my_ndrho = norm_drho(ii)
     310              : 
     311     20673528 :             t7 = 0.1e1_dp/my_rho
     312     20673528 :             t254 = my_ndrho**2
     313     20673528 :             tau_w = t254*t7/0.8e1_dp
     314              : 
     315     20673528 :             IF (my_tau < tau_w) THEN
     316              :                ! enforce z=norm_rho**2/(8._dp*rho*tau) <1
     317      1012581 :                m = 0.5_dp*t254 + 4.0_dp*my_rho*my_tau
     318      1012581 :                my_tau = m/8._dp/my_rho
     319      1012581 :                my_ndrho = SQRT(m)
     320      1012581 :                t254 = m
     321      1012581 :                non_coer = non_coer + 1
     322              :             END IF
     323              : 
     324     20673528 :             t9 = (t6*t7)**(0.1e1_dp/0.3e1_dp)
     325     20673528 :             rs_s1 = t4*t5*t9/0.4e1_dp
     326     20673528 :             phi_s1 = t12/0.2e1_dp
     327     20673528 :             t13 = t1*t12
     328     20673528 :             t14 = pi**2
     329     20673528 :             t16 = (t14*my_rho)**(0.1e1_dp/0.3e1_dp)
     330     20673528 :             k_f_s1 = t13*t16/0.2e1_dp
     331     20673528 :             t19 = SQRT(k_f_s1*t6)
     332     20673528 :             k_s_s1 = 0.2e1_dp*t19
     333     20673528 :             t20 = 0.1e1_dp/phi_s1
     334     20673528 :             t21 = my_ndrho*t20
     335     20673528 :             t22 = 0.1e1_dp/k_s_s1
     336     20673528 :             t_s1 = t21*t22*t7/0.2e1_dp
     337     20673528 :             rs_s2 = rs_s1
     338     20673528 :             phi_s2 = phi_s1
     339     20673528 :             t26 = SQRT(k_f_s1*t6)
     340     20673528 :             k_s_s2 = 0.2e1_dp*t26
     341     20673528 :             t27 = 0.1e1_dp/phi_s2
     342     20673528 :             t28 = my_ndrho*t27
     343     20673528 :             t29 = 0.1e1_dp/k_s_s2
     344     20673528 :             t_s2 = t28*t29*t7/0.2e1_dp
     345     20673528 :             t35 = 0.1e1_dp/A_1
     346     20673528 :             t36 = SQRT(rs_s2)
     347     20673528 :             t39 = t36*rs_s2
     348     20673528 :             t41 = p_1 + 0.1e1_dp
     349     20673528 :             t42 = rs_s2**t41
     350              :             t49 = LOG(0.1e1_dp + t35/(beta_1_1*t36 + beta_2_1*rs_s2 + &
     351              :                                       beta_3_1*t39 + beta_4_1*t42)/0.2e1_dp)
     352     20673528 :             t55 = SQRT(rs_s1)
     353     20673528 :             t58 = t55*rs_s1
     354     20673528 :             t60 = rs_s1**t41
     355              :             t67 = LOG(0.1e1_dp + t35/(beta_1_1*t55 + beta_2_1*rs_s1 + &
     356              :                                       beta_3_1*t58 + beta_4_1*t60)/0.2e1_dp)
     357     20673528 :             t71 = 0.1e1_dp + alpha_1_2*rs_s2
     358     20673528 :             t73 = 0.1e1_dp/A_2
     359     20673528 :             t77 = p_2 + 0.1e1_dp
     360     20673528 :             t78 = rs_s2**t77
     361     20673528 :             t79 = beta_4_2*t78
     362     20673528 :             t80 = beta_1_2*t36 + beta_2_2*rs_s2 + beta_3_2*t39 + t79
     363     20673528 :             t84 = 0.1e1_dp + t73/t80/0.2e1_dp
     364     20673528 :             t85 = LOG(t84)
     365     20673528 :             e_c_u_1_s2 = -0.2e1_dp*A_2*t71*t85
     366     20673528 :             t89 = 0.1e1_dp + alpha_1_2*rs_s1
     367     20673528 :             t94 = rs_s1**t77
     368     20673528 :             t95 = beta_4_2*t94
     369     20673528 :             t96 = beta_1_2*t55 + beta_2_2*rs_s1 + beta_3_2*t58 + t95
     370     20673528 :             t100 = 0.1e1_dp + t73/t96/0.2e1_dp
     371     20673528 :             t101 = LOG(t100)
     372     20673528 :             e_c_u_1_s1 = -0.2e1_dp*A_2*t89*t101
     373     20673528 :             t111 = p_3 + 1._dp
     374     20673528 :             rs = t4*t9/0.4e1_dp
     375     20673528 :             t138 = 0.1e1_dp + alpha_1_1*rs
     376     20673528 :             t140 = SQRT(rs)
     377     20673528 :             t143 = t140*rs
     378     20673528 :             t145 = rs**t41
     379     20673528 :             t146 = beta_4_1*t145
     380     20673528 :             t147 = beta_1_1*t140 + beta_2_1*rs + beta_3_1*t143 + t146
     381     20673528 :             t151 = 0.1e1_dp + t35/t147/0.2e1_dp
     382     20673528 :             t152 = LOG(t151)
     383     20673528 :             e_c_u_0 = -0.2e1_dp*A_1*t138*t152
     384     20673528 :             t161 = rs**t77
     385              :             t168 = LOG(0.1e1_dp + t73/(beta_1_2*t140 + beta_2_2*rs + &
     386              :                                        beta_3_2*t143 + beta_4_2*t161)/0.2e1_dp)
     387     20673528 :             t177 = rs**t111
     388     20673528 :             t186 = 0.1e1_dp/gamma_var
     389     20673528 :             t187 = beta*t186
     390     20673528 :             t189 = phi_s1**2
     391     20673528 :             t190 = t189*phi_s1
     392     20673528 :             t191 = 0.1e1_dp/t190
     393     20673528 :             t193 = EXP(-e_c_u_1_s1*t186*t191)
     394     20673528 :             t194 = t193 - 0.1e1_dp
     395     20673528 :             A_s1 = t187/t194
     396     20673528 :             t196 = gamma_var*t190
     397     20673528 :             t197 = t_s1**2
     398     20673528 :             t198 = A_s1*t197
     399     20673528 :             t199 = 0.1e1_dp + t198
     400     20673528 :             t201 = A_s1**2
     401     20673528 :             t202 = t197**2
     402     20673528 :             t204 = 0.1e1_dp + t198 + t201*t202
     403     20673528 :             t205 = 0.1e1_dp/t204
     404     20673528 :             t208 = 0.1e1_dp + t187*t197*t199*t205
     405     20673528 :             t209 = LOG(t208)
     406     20673528 :             epsilon_cGGA_1_0 = e_c_u_1_s1 + t196*t209
     407     20673528 :             t211 = phi_s2**2
     408     20673528 :             t212 = t211*phi_s2
     409     20673528 :             t213 = 0.1e1_dp/t212
     410     20673528 :             t215 = EXP(-e_c_u_1_s2*t186*t213)
     411     20673528 :             t216 = t215 - 0.1e1_dp
     412     20673528 :             A_s2 = t187/t216
     413     20673528 :             t218 = gamma_var*t212
     414     20673528 :             t219 = t_s2**2
     415     20673528 :             t220 = A_s2*t219
     416     20673528 :             t221 = t220 + 0.1e1_dp
     417     20673528 :             t223 = A_s2**2
     418     20673528 :             t224 = t219**2
     419     20673528 :             t226 = 0.1e1_dp + t220 + t223*t224
     420     20673528 :             t227 = 0.1e1_dp/t226
     421     20673528 :             t230 = 0.1e1_dp + t187*t219*t221*t227
     422     20673528 :             t231 = LOG(t230)
     423     20673528 :             epsilon_cGGA_0_1 = e_c_u_1_s2 + t218*t231
     424     20673528 :             t233 = SQRT(t1*t16*t6)
     425     20673528 :             k_s = 0.2e1_dp*t233
     426     20673528 :             t234 = 0.1e1_dp/k_s
     427     20673528 :             t235 = my_ndrho*t234
     428     20673528 :             t = t235*t7/0.2e1_dp
     429     20673528 :             t238 = EXP(-e_c_u_0*t186)
     430     20673528 :             t239 = -0.1e1_dp + t238
     431     20673528 :             A = t187/t239
     432     20673528 :             t241 = t**2
     433     20673528 :             t242 = A*t241
     434     20673528 :             t243 = 0.1e1_dp + t242
     435     20673528 :             t245 = A**2
     436     20673528 :             t246 = t241**2
     437     20673528 :             t248 = 0.1e1_dp + t242 + t245*t246
     438     20673528 :             t249 = 0.1e1_dp/t248
     439     20673528 :             t252 = 0.1e1_dp + t187*t241*t243*t249
     440     20673528 :             t253 = LOG(t252)
     441     20673528 :             epsilon_cGGA = e_c_u_0 + gamma_var*t253
     442     20673528 :             ma = MAX(epsilon_cGGA_1_0, epsilon_cGGA)
     443     20673528 :             mb = MAX(epsilon_cGGA_0_1, epsilon_cGGA)
     444     20673528 :             t256 = tau_w**2
     445     20673528 :             t260 = ma/0.2e1_dp + mb/0.2e1_dp
     446     20673528 :             t263 = 0.53e0_dp*epsilon_cGGA*t256 - 0.153e1_dp*t256*t260
     447     20673528 :             t264 = my_tau**2
     448     20673528 :             t265 = 0.1e1_dp/t264
     449     20673528 :             epsilon_cRevPKZB = epsilon_cGGA + t263*t265
     450     20673528 :             t267 = my_rho*epsilon_cRevPKZB
     451     20673528 :             t268 = d*epsilon_cRevPKZB
     452     20673528 :             t269 = t256*tau_w
     453     20673528 :             t271 = 0.1e1_dp/t264/my_tau
     454     20673528 :             t272 = t269*t271
     455     20673528 :             t274 = 0.1e1_dp + t268*t272
     456     20673528 :             t275 = t254*t1
     457     20673528 :             t276 = t14**(0.1e1_dp/0.3e1_dp)
     458     20673528 :             t277 = t276**2
     459     20673528 :             t278 = 0.1e1_dp/t277
     460     20673528 :             t279 = my_rho**2
     461     20673528 :             t280 = my_rho**(0.1e1_dp/0.3e1_dp)
     462     20673528 :             t281 = t280**2
     463     20673528 :             t284 = t278/t281/t279
     464     20673528 :             p = t275*t284/0.12e2_dp
     465     20673528 :             t286 = 0.1e1_dp/my_tau
     466     20673528 :             z = tau_w*t286
     467     20673528 :             t288 = 0.1e1_dp/z - 0.1e1_dp
     468     20673528 :             alpha = 0.5e1_dp/0.3e1_dp*p*t288
     469     20673528 :             t290 = alpha - 0.1e1_dp
     470     20673528 :             t291 = b*alpha
     471     20673528 :             t293 = 0.1e1_dp + t291*t290
     472     20673528 :             t294 = SQRT(t293)
     473     20673528 :             t295 = 0.1e1_dp/t294
     474     20673528 :             tildeq_b = 0.9e1_dp/0.20e2_dp*t290*t295 + 0.2e1_dp/0.3e1_dp*p
     475     20673528 :             t299 = z**2
     476     20673528 :             t301 = 0.1e1_dp + t299
     477     20673528 :             t302 = t301**2
     478     20673528 :             t303 = 0.1e1_dp/t302
     479     20673528 :             t305 = 0.10e2_dp/0.81e2_dp + c*t299*t303
     480     20673528 :             t307 = tildeq_b**2
     481     20673528 :             t310 = p**2
     482     20673528 :             t313 = SQRT(0.18e2_dp*t299 + 0.50e2_dp*t310)
     483     20673528 :             t316 = 0.1e1_dp/kappa
     484     20673528 :             t319 = SQRT(e_var)
     485     20673528 :             t322 = e_var*mu
     486              :             t325 = t305*p + 0.146e3_dp/0.2025e4_dp*t307 - 0.73e2_dp/ &
     487              :                    0.4050e4_dp*tildeq_b*t313 + 0.100e3_dp/0.6561e4_dp*t316* &
     488     20673528 :                    t310 + 0.4e1_dp/0.45e2_dp*t319*t299 + t322*t310*p
     489     20673528 :             t327 = 0.1e1_dp + t319*p
     490     20673528 :             t328 = t327**2
     491     20673528 :             t329 = 0.1e1_dp/t328
     492     20673528 :             t331 = 0.1e1_dp + t325*t329*t316
     493     20673528 :             Fx = 0.1e1_dp + kappa - kappa/t331
     494     20673528 :             ex_unif = -0.3e1_dp/0.4e1_dp*t1*t16*t6
     495     20673528 :             t337 = my_rho*ex_unif
     496              : 
     497     20673528 :             IF (grad_deriv >= 0) THEN
     498              :                e_0(ii) = e_0(ii) + &
     499     20673528 :                          scale_ec*t267*t274 + scale_ex*t337*Fx
     500              :             END IF
     501              : 
     502     20673528 :             IF (abs_grad_deriv > 0) THEN
     503     20673528 :                t340 = t9**2
     504     20673528 :                t343 = 0.1e1_dp/t279
     505     20673528 :                t344 = 0.1e1_dp/t340*t6*t343
     506     20673528 :                rsrho = -t4*t344/0.12e2_dp
     507     20673528 :                t351 = t147**2
     508              :                e_c_u_0rho = -0.2e1_dp*A_1*alpha_1_1*rsrho*t152 + t138/ &
     509              :                             t351*(beta_1_1/t140*rsrho/0.2e1_dp + beta_2_1*rsrho + &
     510              :                                   0.3e1_dp/0.2e1_dp*beta_3_1*t140*rsrho + t146*t41*rsrho/ &
     511     20673528 :                                   rs)/t151
     512     20673528 :                t370 = t16**2
     513     20673528 :                t371 = 0.1e1_dp/t370
     514     20673528 :                t376 = k_s**2
     515              :                trho = -my_ndrho/t376*t7/t233*t1*t371*t14*t6 &
     516     20673528 :                       /0.6e1_dp - t235*t343/0.2e1_dp
     517     20673528 :                t383 = gamma_var**2
     518     20673528 :                t385 = beta/t383
     519     20673528 :                t386 = t239**2
     520     20673528 :                Arho = t385/t386*e_c_u_0rho*t238
     521     20673528 :                t390 = t187*t
     522     20673528 :                t391 = t243*t249
     523     20673528 :                t395 = Arho*t241
     524     20673528 :                t396 = A*t
     525     20673528 :                t398 = 0.2e1_dp*t396*trho
     526     20673528 :                t403 = t187*t241
     527     20673528 :                t404 = t248**2
     528     20673528 :                t406 = t243/t404
     529     20673528 :                t410 = t241*t
     530     20673528 :                t411 = t245*t410
     531     20673528 :                t419 = 0.1e1_dp/t252
     532              :                epsilon_cGGArho = e_c_u_0rho + gamma_var*(0.2e1_dp*t390*t391 &
     533              :                                                          *trho + t187*t241*(t395 + t398)*t249 - t403*t406*(t395 + &
     534     20673528 :                                                                              t398 + 0.2e1_dp*A*t246*Arho + 0.4e1_dp*t411*trho))*t419
     535     20673528 :                tau_wrho = -t254*t343/0.8e1_dp
     536     20673528 :                prho = -0.2e1_dp/0.9e1_dp*t275*t278/t281/t279/my_rho
     537     20673528 :                zrho = tau_wrho*t286
     538     20673528 :                t430 = p/t299
     539              :                alpharho = 0.5e1_dp/0.3e1_dp*prho*t288 - 0.5e1_dp/0.3e1_dp &
     540     20673528 :                           *t430*zrho
     541     20673528 :                t437 = t290/t294/t293
     542              :                tildeq_brho = 0.9e1_dp/0.20e2_dp*alpharho*t295 - 0.9e1_dp/ &
     543              :                              0.40e2_dp*t437*(b*alpharho*t290 + t291*alpharho) + &
     544     20673528 :                              0.2e1_dp/0.3e1_dp*prho
     545     20673528 :                t445 = c*z
     546     20673528 :                t450 = c*t299*z
     547     20673528 :                t452 = 0.1e1_dp/t302/t301
     548     20673528 :                t464 = tildeq_b/t313
     549     20673528 :                t472 = t316*p
     550     20673528 :                t475 = t319*z
     551     20673528 :                t485 = t325/t328/t327
     552     20673528 :                t489 = t331**2
     553     20673528 :                t490 = 0.1e1_dp/t489
     554     20673528 :                rs_s1rho = -t4*t5*t344/0.12e2_dp
     555     20673528 :                k_f_s1rho = t13*t371*t14/0.6e1_dp
     556     20673528 :                t505 = k_s_s1**2
     557              :                t_s1rho = -t21/t505*t7/t19*k_f_s1rho*t6/0.2e1_dp - t21 &
     558     20673528 :                          *t22*t343/0.2e1_dp
     559     20673528 :                t513 = A_2*alpha_1_2
     560     20673528 :                t517 = t96**2
     561              :                e_c_u_1_s1rho = -0.2e1_dp*t513*rs_s1rho*t101 + t89/t517* &
     562              :                                (beta_1_2/t55*rs_s1rho/0.2e1_dp + beta_2_2*rs_s1rho + &
     563              :                                 0.3e1_dp/0.2e1_dp*beta_3_2*t55*rs_s1rho + t95*t77* &
     564     20673528 :                                 rs_s1rho/rs_s1)/t100
     565     20673528 :                t536 = t194**2
     566     20673528 :                A_s1rho = t385/t536*e_c_u_1_s1rho*t191*t193
     567     20673528 :                t541 = t187*t_s1
     568     20673528 :                t542 = t199*t205
     569     20673528 :                t546 = A_s1rho*t197
     570     20673528 :                t547 = A_s1*t_s1
     571     20673528 :                t549 = 0.2e1_dp*t547*t_s1rho
     572     20673528 :                t554 = t187*t197
     573     20673528 :                t555 = t204**2
     574     20673528 :                t557 = t199/t555
     575     20673528 :                t561 = t197*t_s1
     576     20673528 :                t562 = t201*t561
     577     20673528 :                t569 = 0.1e1_dp/t208
     578     20673528 :                t571 = epsilon_cGGA .LT. epsilon_cGGA_1_0
     579     20673528 :                IF (t571) THEN
     580              :                   marho = e_c_u_1_s1rho + t196*(0.2e1_dp*t541*t542 &
     581              :                                                 *t_s1rho + t187*t197*(t546 + t549)*t205 - t554*t557*(t546 &
     582              :                                                                               + t549 + 0.2e1_dp*A_s1*t202*A_s1rho + 0.4e1_dp*t562* &
     583     20619872 :                                                                                                      t_s1rho))*t569
     584              :                ELSE
     585              :                   marho = epsilon_cGGArho
     586              :                END IF
     587     20673528 :                rs_s2rho = rs_s1rho
     588     20673528 :                t574 = k_s_s2**2
     589              :                t_s2rho = -t28/t574*t7/t26*k_f_s1rho*t6/0.2e1_dp - t28 &
     590     20673528 :                          *t29*t343/0.2e1_dp
     591     20673528 :                t585 = t80**2
     592              :                e_c_u_1_s2rho = -0.2e1_dp*t513*rs_s2rho*t85 + t71/t585*( &
     593              :                                beta_1_2/t36*rs_s2rho/0.2e1_dp + beta_2_2*rs_s2rho + &
     594              :                                0.3e1_dp/0.2e1_dp*beta_3_2*t36*rs_s2rho + t79*t77* &
     595     20673528 :                                rs_s2rho/rs_s2)/t84
     596     20673528 :                t604 = t216**2
     597     20673528 :                A_s2rho = t385/t604*e_c_u_1_s2rho*t213*t215
     598     20673528 :                t609 = t187*t_s2
     599     20673528 :                t610 = t221*t227
     600     20673528 :                t614 = A_s2rho*t219
     601     20673528 :                t615 = A_s2*t_s2
     602     20673528 :                t617 = 0.2e1_dp*t615*t_s2rho
     603     20673528 :                t622 = t187*t219
     604     20673528 :                t623 = t226**2
     605     20673528 :                t625 = t221/t623
     606     20673528 :                t629 = t219*t_s2
     607     20673528 :                t630 = t223*t629
     608     20673528 :                t637 = 0.1e1_dp/t230
     609     20673528 :                t639 = epsilon_cGGA .LT. epsilon_cGGA_0_1
     610              :                IF (t639) THEN
     611              :                   mbrho = e_c_u_1_s2rho + t218*(0.2e1_dp*t609*t610 &
     612              :                                                 *t_s2rho + t187*t219*(t614 + t617)*t227 - t622*t625*(t614 &
     613              :                                                                               + t617 + 0.2e1_dp*A_s2*t224*A_s2rho + 0.4e1_dp*t630* &
     614     20673528 :                                                                                                      t_s2rho))*t637
     615              :                ELSE
     616              :                   mbrho = epsilon_cGGArho
     617              :                END IF
     618     20673528 :                t642 = epsilon_cGGA*tau_w
     619     20673528 :                t645 = tau_w*t260
     620              :                epsilon_cRevPKZBrho = epsilon_cGGArho + (0.53e0_dp* &
     621              :                                                         epsilon_cGGArho*t256 + 0.106e1_dp*t642*tau_wrho - 0.306e1_dp* &
     622              :                                                         t645*tau_wrho - 0.153e1_dp*t256*(marho/0.2e1_dp + mbrho/ &
     623     20673528 :                                                                                          0.2e1_dp))*t265
     624     20673528 :                t659 = t256*t271
     625              : 
     626     20673528 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     627              :                   e_rho(ii) = e_rho(ii) + &
     628              :                               scale_ec*(epsilon_cRevPKZB*t274 + my_rho* &
     629              :                                         epsilon_cRevPKZBrho*t274 + t267*(d*epsilon_cRevPKZBrho*t272 &
     630              :                                                                          + 0.3e1_dp*t268*t659*tau_wrho)) + scale_ex*(ex_unif*Fx - &
     631              :                                                                                              my_rho*pi*t1*t371*Fx/0.4e1_dp + t337* &
     632              :                                                                              t490*(((0.2e1_dp*t445*t303*zrho - 0.4e1_dp*t450*t452* &
     633              :                                                                             zrho)*p + t305*prho + 0.292e3_dp/0.2025e4_dp*tildeq_b* &
     634              :                                                                             tildeq_brho - 0.73e2_dp/0.4050e4_dp*tildeq_brho*t313 - &
     635              :                                                                          0.73e2_dp/0.8100e4_dp*t464*(0.36e2_dp*z*zrho + 0.100e3_dp &
     636              :                                                                            *p*prho) + 0.200e3_dp/0.6561e4_dp*t472*prho + 0.8e1_dp/ &
     637              :                                                                              0.45e2_dp*t475*zrho + 0.3e1_dp*t322*t310*prho)*t329 - &
     638     20673528 :                                                                                                            0.2e1_dp*t485*t319*prho))
     639              :                END IF
     640              : 
     641     20673528 :                tnorm_drho = t234*t7/0.2e1_dp
     642              :                Hnorm_drho = gamma_var*(0.2e1_dp*t390*t391*tnorm_drho + &
     643              :                                        0.2e1_dp*t187*t410*A*tnorm_drho*t249 - t403*t406*( &
     644     20673528 :                                        0.2e1_dp*t396*tnorm_drho + 0.4e1_dp*t411*tnorm_drho))*t419
     645     20673528 :                tau_wnorm_drho = my_ndrho*t7/0.4e1_dp
     646     20673528 :                pnorm_drho = my_ndrho*t1*t284/0.6e1_dp
     647     20673528 :                znorm_drho = tau_wnorm_drho*t286
     648              :                alphanorm_drho = 0.5e1_dp/0.3e1_dp*pnorm_drho*t288 - &
     649     20673528 :                                 0.5e1_dp/0.3e1_dp*t430*znorm_drho
     650              :                tildeq_bnorm_drho = 0.9e1_dp/0.20e2_dp*alphanorm_drho*t295 - &
     651              :                                    0.9e1_dp/0.40e2_dp*t437*(b*alphanorm_drho*t290 + t291* &
     652     20673528 :                                                             alphanorm_drho) + 0.2e1_dp/0.3e1_dp*pnorm_drho
     653     20673528 :                t_s1norm_drho = t20*t22*t7/0.2e1_dp
     654     20673528 :                IF (t571) THEN
     655              :                   manorm_drho = t196*(0.2e1_dp*t541*t542* &
     656              :                                       t_s1norm_drho + 0.2e1_dp*t187*t561*A_s1*t_s1norm_drho*t205 &
     657              :                                       - t554*t557*(0.2e1_dp*t547*t_s1norm_drho + 0.4e1_dp*t562 &
     658     20619872 :                                                    *t_s1norm_drho))*t569
     659              :                ELSE
     660              :                   manorm_drho = Hnorm_drho
     661              :                END IF
     662     20673528 :                t_s2norm_drho = t27*t29*t7/0.2e1_dp
     663              :                IF (t639) THEN
     664              :                   mbnorm_drho = t218*(0.2e1_dp*t609*t610* &
     665              :                                       t_s2norm_drho + 0.2e1_dp*t187*t629*A_s2*t_s2norm_drho*t227 &
     666              :                                       - t622*t625*(0.2e1_dp*t615*t_s2norm_drho + 0.4e1_dp*t630 &
     667     20673528 :                                                    *t_s2norm_drho))*t637
     668              :                ELSE
     669              :                   mbnorm_drho = Hnorm_drho
     670              :                END IF
     671              :                epsilon_cRevPKZBnorm_drho = Hnorm_drho + (0.53e0_dp*Hnorm_drho* &
     672              :                                                          t256 + 0.106e1_dp*t642*tau_wnorm_drho - 0.306e1_dp*t645* &
     673              :                                                          tau_wnorm_drho - 0.153e1_dp*t256*(manorm_drho/0.2e1_dp + &
     674     20673528 :                                                                                            mbnorm_drho/0.2e1_dp))*t265
     675              : 
     676     20673528 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     677              :                   e_ndrho(ii) = e_ndrho(ii) + &
     678              :                                 scale_ec*(my_rho*epsilon_cRevPKZBnorm_drho* &
     679              :                                           t274 + t267*(d*epsilon_cRevPKZBnorm_drho*t272 + 0.3e1_dp* &
     680              :                                                        t268*t659*tau_wnorm_drho)) + scale_ex*t337*t490*((( &
     681              :                                                                                0.2e1_dp*t445*t303*znorm_drho - 0.4e1_dp*t450*t452* &
     682              :                                                                          znorm_drho)*p + t305*pnorm_drho + 0.292e3_dp/0.2025e4_dp* &
     683              :                                                                                tildeq_b*tildeq_bnorm_drho - 0.73e2_dp/0.4050e4_dp* &
     684              :                                                                              tildeq_bnorm_drho*t313 - 0.73e2_dp/0.8100e4_dp*t464*( &
     685              :                                                                                0.36e2_dp*z*znorm_drho + 0.100e3_dp*p*pnorm_drho) + &
     686              :                                                                        0.200e3_dp/0.6561e4_dp*t472*pnorm_drho + 0.8e1_dp/0.45e2_dp &
     687              :                                                                           *t475*znorm_drho + 0.3e1_dp*t322*t310*pnorm_drho)*t329 - &
     688     20673528 :                                                                                                       0.2e1_dp*t485*t319*pnorm_drho)
     689              :                END IF
     690              : 
     691     20673528 :                epsilon_cRevPKZBtau = -0.2e1_dp*t263*t271
     692     20673528 :                t799 = t264**2
     693     20673528 :                ztau = -tau_w*t265
     694     20673528 :                alphatau = -0.5e1_dp/0.3e1_dp*t430*ztau
     695              :                tildeq_btau = 0.9e1_dp/0.20e2_dp*alphatau*t295 - 0.9e1_dp/ &
     696     20673528 :                              0.40e2_dp*t437*(b*alphatau*t290 + t291*alphatau)
     697              : 
     698     20673528 :                IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
     699              :                   e_tau(ii) = e_tau(ii) + &
     700              :                               scale_ec*(my_rho*epsilon_cRevPKZBtau*t274 + t267* &
     701              :                                         (d*epsilon_cRevPKZBtau*t272 - 0.3e1_dp*t268*t269/t799)) + &
     702              :                               scale_ex*t337*t490*((0.2e1_dp*t445*t303*ztau - 0.4e1_dp &
     703              :                                                    *t450*t452*ztau)*p + 0.292e3_dp/0.2025e4_dp*tildeq_b* &
     704              :                                                   tildeq_btau - 0.73e2_dp/0.4050e4_dp*tildeq_btau*t313 - &
     705              :                                                   0.73e2_dp/0.225e3_dp*t464*z*ztau + 0.8e1_dp/0.45e2_dp* &
     706     20673528 :                                                   t475*ztau)*t329
     707              :                END IF
     708              :             END IF
     709              :          END IF
     710              :       END DO
     711              : 
     712              : !$OMP     END DO
     713              : 
     714              :    END SUBROUTINE tpss_lda_calc
     715              : 
     716              : END MODULE xc_tpss
        

Generated by: LCOV version 2.0-1