3 ! { dg-options "-O1 -ffast-math" }
6 INTEGER, PARAMETER :: dp
=8
8 PUBLIC
:: b97_lda_info
, b97_lsd_info
, b97_lda_eval
, b97_lsd_eval
10 SUBROUTINE b97_lsd_eval(rho_set
,deriv_set
,grad_deriv
,b97_params
)
11 INTEGER, INTENT(in
) :: grad_deriv
12 INTEGER :: handle
, npoints
, param
, stat
14 REAL(kind
=dp
) :: epsilon_drho
, epsilon_rho
, &
16 REAL(kind
=dp
), DIMENSION(:, :, :), POINTER :: dummy
, e_0
, e_ndra
, &
17 e_ndra_ndra
, e_ndra_ndrb
, e_ndra_ra
, e_ndra_rb
, e_ndrb
, e_ndrb_ndrb
, &
18 e_ndrb_ra
, e_ndrb_rb
, e_ra
, e_ra_ra
, e_ra_rb
, e_rb
, e_rb_rb
, &
19 norm_drhoa
, norm_drhob
, rhoa
, rhob
20 IF (.NOT
. failure
) THEN
22 rhoa
=rhoa
, rhob
=rhob
, norm_drhoa
=norm_drhoa
,&
23 norm_drhob
=norm_drhob
, e_0
=e_0
, &
24 e_ra
=e_ra
, e_rb
=e_rb
, &
25 e_ndra
=e_ndra
, e_ndrb
=e_ndrb
, &
26 e_ra_ra
=e_ra_ra
, e_ra_rb
=e_ra_rb
, e_rb_rb
=e_rb_rb
,&
27 e_ra_ndra
=e_ndra_ra
, e_ra_ndrb
=e_ndrb_ra
, &
28 e_rb_ndrb
=e_ndrb_rb
, e_rb_ndra
=e_ndra_rb
,&
29 e_ndra_ndra
=e_ndra_ndra
, e_ndrb_ndrb
=e_ndrb_ndrb
,&
30 e_ndra_ndrb
=e_ndra_ndrb
,&
31 grad_deriv
=grad_deriv
, npoints
=npoints
, &
32 epsilon_rho
=epsilon_rho
,epsilon_drho
=epsilon_drho
,&
33 param
=param
,scale_c_in
=scale_c
,scale_x_in
=scale_x
)
35 END SUBROUTINE b97_lsd_eval
36 SUBROUTINE b97_lsd_calc(rhoa
, rhob
, norm_drhoa
, norm_drhob
,&
37 e_0
, e_ra
, e_rb
, e_ndra
, e_ndrb
, &
38 e_ra_ndra
,e_ra_ndrb
, e_rb_ndra
, e_rb_ndrb
,&
39 e_ndra_ndra
, e_ndrb_ndrb
, e_ndra_ndrb
, &
40 e_ra_ra
, e_ra_rb
, e_rb_rb
,&
41 grad_deriv
,npoints
,epsilon_rho
,epsilon_drho
, &
42 param
, scale_c_in
, scale_x_in
)
43 REAL(kind
=dp
), DIMENSION(*), INTENT(in
) :: rhoa
, rhob
, norm_drhoa
, &
45 REAL(kind
=dp
), DIMENSION(*), INTENT(inout
) :: e_0
, e_ra
, e_rb
, e_ndra
, &
46 e_ndrb
, e_ra_ndra
, e_ra_ndrb
, e_rb_ndra
, e_rb_ndrb
, e_ndra_ndra
, &
47 e_ndrb_ndrb
, e_ndra_ndrb
, e_ra_ra
, e_ra_rb
, e_rb_rb
48 INTEGER, INTENT(in
) :: grad_deriv
, npoints
49 REAL(kind
=dp
), INTENT(in
) :: epsilon_rho
, epsilon_drho
50 INTEGER, INTENT(in
) :: param
51 REAL(kind
=dp
), INTENT(in
) :: scale_c_in
, scale_x_in
52 REAL(kind
=dp
) :: A_1
, A_2
, A_3
, alpha_1_1
, alpha_1_2
, alpha_1_3
, alpha_c
, &
53 t133
, t134
, t1341
, t1348
, t1351
, t1360
, t1368
, t138
, t1388
, t139
, &
54 u_x_bnorm_drhobnorm_drhob
, u_x_brhob
, u_x_brhobnorm_drhob
, u_x_brhobrhob
55 SELECT
CASE(grad_deriv
)
58 IF (rho
>epsilon_rho
) THEN
59 IF (grad_deriv
/=0) THEN
60 IF (grad_deriv
>1 .OR
. grad_deriv
<-1) THEN
61 alpha_c1rhob
= alpha_crhob
63 t1360
= -0.4e1_dp
* t105
* t290
* chirhobrhob
+ (-0.2e1_dp
* t239
&
64 * t257
+ t709
* t1236
* t711
* t62
/ 0.2e1_dp
- e_c_u_0rhobrhob
) * f
&
65 * t108
+ t438
* f1rhob
* t108
+ 0.4e1_dp
* t439
* t443
+ t1341
* &
66 0.4e1_dp
* t1348
* t443
+ 0.4e1_dp
* t1351
* t443
+ 0.12e2_dp
* t113
&
67 * t107
* t1299
+ 0.4e1_dp
* t113
* t289
* chirhobrhob
68 IF (grad_deriv
>1 .OR
. grad_deriv
==-2) THEN
69 exc_rhob_rhob
= scale_x
* (-t4
* t6
/ t1152
* gx_b
/ &
70 0.6e1_dp
+ e_lsda_x_brhob
* (u_x_b1rhob
* t31
+ u_x_b
* u_x_b1rhob
*&
71 u_x_brhobrhob
* c_x_2
)) + scale_c
* (((e_c_u_0rhobrhob
+ (0.2e1_dp
*&
72 t726
* t1270
* t278
- t266
* (-t731
* t1205
/ 0.4e1_dp
+ t267
* &
73 t1205
* t647
) * t278
- t757
* t1270
* t759
* t80
/ 0.2e1_dp
) * f
* &
74 t110
+ alpha_crhob
* f1rhob
* t110
- 0.4e1_dp
* t431
* t435
+ &
75 alpha_c1rhob
* frhob
* t110
+ alpha_c
* frhobrhob
* t110
- 0.4e1_dp
&
76 * t433
* t435
- 0.4e1_dp
* t1321
* t435
- 0.4e1_dp
* t1324
* t435
- &
77 0.12e2_dp
* t105
* t796
* t1299
+ t1360
) * rho
+ epsilon_c_unifrhob
&
79 e_rb_rb(ii
)=e_rb_rb(ii
)+exc_rhob_rhob
83 END IF ! rho>epsilon_rho
86 END SUBROUTINE b97_lsd_calc