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 calculates fxc in the spirit of the b97 exchange/correlation functional
      10             : !> \author jgh
      11             : ! **************************************************************************************************
      12             : MODULE xc_b97_fxc
      13             :    USE kinds,                           ONLY: dp
      14             :    USE pw_types,                        ONLY: pw_r3d_rs_type
      15             : #include "../base/base_uses.f90"
      16             : 
      17             :    IMPLICIT NONE
      18             :    PRIVATE
      19             : 
      20             :    PUBLIC :: b97_fxc_eval, b97_fcc_eval
      21             : 
      22             : CONTAINS
      23             : 
      24             : ! **************************************************************************************************
      25             : !> \brief ...
      26             : !> \param rhos ...
      27             : !> \param norm_drhos ...
      28             : !> \param fxc ...
      29             : !> \param gx ...
      30             : !> \param cx ...
      31             : !> \param eps_rho ...
      32             : ! **************************************************************************************************
      33          16 :    SUBROUTINE b97_fxc_eval(rhos, norm_drhos, fxc, gx, cx, eps_rho)
      34             :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: rhos, norm_drhos
      35             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: fxc
      36             :       REAL(KIND=dp), INTENT(IN)                          :: gx
      37             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: cx
      38             :       REAL(KIND=dp), INTENT(IN)                          :: eps_rho
      39             : 
      40             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'b97_fxc_eval'
      41             : 
      42             :       INTEGER                                            :: handle, i, io, j, k, norder
      43             :       INTEGER, DIMENSION(2, 3)                           :: bo
      44             :       REAL(KIND=dp)                                      :: drho, gval, rho, sval, ux
      45             : 
      46          16 :       CALL timeset(routineN, handle)
      47             : 
      48          16 :       norder = SIZE(cx)
      49         160 :       bo(1:2, 1:3) = rhos%pw_grid%bounds_local(1:2, 1:3)
      50             : !$OMP PARALLEL DO PRIVATE(i,j,k,io,rho,drho,sval,gval,ux) DEFAULT(NONE)&
      51          16 : !$OMP SHARED(bo,rhos,norm_drhos,fxc,gx,cx,eps_rho,norder)
      52             :       DO k = bo(1, 3), bo(2, 3)
      53             :          DO j = bo(1, 2), bo(2, 2)
      54             :             DO i = bo(1, 1), bo(2, 1)
      55             : 
      56             :                rho = rhos%array(i, j, k)
      57             :                drho = norm_drhos%array(i, j, k)
      58             :                IF (rho > eps_rho) THEN
      59             :                   sval = gx*(drho/rho**1.33333333333333_dp)**2
      60             :                   ux = sval/(1._dp + sval)
      61             :                   gval = 0.0_dp
      62             :                   DO io = 0, norder - 1
      63             :                      gval = gval + cx(io + 1)*(ux**io)
      64             :                   END DO
      65             :                   fxc%array(i, j, k) = fxc%array(i, j, k)*gval
      66             :                END IF
      67             : 
      68             :             END DO
      69             :          END DO
      70             :       END DO
      71             : 
      72          16 :       CALL timestop(handle)
      73             : 
      74          16 :    END SUBROUTINE b97_fxc_eval
      75             : 
      76             : ! **************************************************************************************************
      77             : !> \brief ...
      78             : !> \param rhoa ...
      79             : !> \param rhob ...
      80             : !> \param norm_drhoa ...
      81             : !> \param norm_drhob ...
      82             : !> \param fcc ...
      83             : !> \param gcc ...
      84             : !> \param cco ...
      85             : !> \param eps_rho ...
      86             : ! **************************************************************************************************
      87           4 :    SUBROUTINE b97_fcc_eval(rhoa, rhob, norm_drhoa, norm_drhob, fcc, gcc, cco, eps_rho)
      88             :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: rhoa, rhob, norm_drhoa, norm_drhob
      89             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: fcc
      90             :       REAL(KIND=dp), INTENT(IN)                          :: gcc
      91             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: cco
      92             :       REAL(KIND=dp), INTENT(IN)                          :: eps_rho
      93             : 
      94             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'b97_fcc_eval'
      95             : 
      96             :       INTEGER                                            :: handle, i, io, j, k, norder
      97             :       INTEGER, DIMENSION(2, 3)                           :: bo
      98             :       REAL(KIND=dp)                                      :: dra, drb, gval, ra, rb, sa, sb, sval, ux
      99             : 
     100           4 :       CALL timeset(routineN, handle)
     101             : 
     102           4 :       norder = SIZE(cco)
     103          40 :       bo(1:2, 1:3) = rhoa%pw_grid%bounds_local(1:2, 1:3)
     104             : !$OMP PARALLEL DO PRIVATE(i,j,k,ra,rb,dra,drb,sa,sb,sval,gval,ux,io) DEFAULT(NONE)&
     105           4 : !$OMP SHARED(bo,rhoa,rhob,norm_drhoa,norm_drhob,fcc,gcc,cco,norder,eps_rho)
     106             :       DO k = bo(1, 3), bo(2, 3)
     107             :          DO j = bo(1, 2), bo(2, 2)
     108             :             DO i = bo(1, 1), bo(2, 1)
     109             : 
     110             :                ra = rhoa%array(i, j, k)
     111             :                rb = rhob%array(i, j, k)
     112             :                dra = norm_drhoa%array(i, j, k)
     113             :                drb = norm_drhob%array(i, j, k)
     114             :                IF (ra > eps_rho .AND. rb > eps_rho) THEN
     115             :                   sa = (dra/ra**1.33333333333333_dp)**2
     116             :                   sb = (drb/rb**1.33333333333333_dp)**2
     117             :                   sval = 0.5_dp*gcc*(sa + sb)
     118             :                   ux = sval/(1._dp + sval)
     119             :                   gval = 0.0_dp
     120             :                   DO io = 0, norder - 1
     121             :                      gval = gval + cco(io + 1)*(ux**io)
     122             :                   END DO
     123             :                   fcc%array(i, j, k) = fcc%array(i, j, k)*gval
     124             :                END IF
     125             : 
     126             :             END DO
     127             :          END DO
     128             :       END DO
     129             : 
     130           4 :       CALL timestop(handle)
     131             : 
     132           4 :    END SUBROUTINE b97_fcc_eval
     133             : 
     134             : END MODULE xc_b97_fxc