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 the lyp correlation functional
10 : !> \par History
11 : !> 11.2003 created [fawzi]
12 : !> \author fawzi
13 : ! **************************************************************************************************
14 : MODULE xc_lyp
15 : USE bibliography, ONLY: Lee1988,&
16 : cite_reference
17 : USE input_section_types, ONLY: section_vals_type,&
18 : section_vals_val_get
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: pi
21 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
22 : deriv_norm_drhoa,&
23 : deriv_norm_drhob,&
24 : deriv_rho,&
25 : deriv_rhoa,&
26 : deriv_rhob
27 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
28 : xc_dset_get_derivative
29 : USE xc_derivative_types, ONLY: xc_derivative_get,&
30 : xc_derivative_type
31 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
32 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
33 : xc_rho_set_type
34 : #include "../base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 : PRIVATE
38 :
39 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_lyp'
41 : REAL(kind=dp), PARAMETER, PRIVATE :: a = 0.04918_dp, b = 0.132_dp, &
42 : c = 0.2533_dp, d = 0.349_dp
43 :
44 : PUBLIC :: lyp_lda_info, lyp_lsd_info, lyp_lda_eval, lyp_lsd_eval
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief return various information on the functional
49 : !> \param reference string with the reference of the actual functional
50 : !> \param shortform string with the shortform of the functional name
51 : !> \param needs the components needed by this functional are set to
52 : !> true (does not set the unneeded components to false)
53 : !> \param max_deriv ...
54 : !> \par History
55 : !> 11.2003 created [fawzi]
56 : !> \author fawzi
57 : ! **************************************************************************************************
58 7608 : SUBROUTINE lyp_lda_info(reference, shortform, needs, max_deriv)
59 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
60 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
61 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
62 :
63 7608 : IF (PRESENT(reference)) THEN
64 87 : reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}"
65 : END IF
66 7608 : IF (PRESENT(shortform)) THEN
67 87 : shortform = "Lee-Yang-Parr correlation energy functional (LDA)"
68 : END IF
69 7608 : IF (PRESENT(needs)) THEN
70 7521 : needs%rho = .TRUE.
71 7521 : needs%rho_1_3 = .TRUE.
72 7521 : needs%norm_drho = .TRUE.
73 : END IF
74 7608 : IF (PRESENT(max_deriv)) max_deriv = 3
75 :
76 7608 : END SUBROUTINE lyp_lda_info
77 :
78 : ! **************************************************************************************************
79 : !> \brief return various information on the functional
80 : !> \param reference string with the reference of the actual functional
81 : !> \param shortform string with the shortform of the functional name
82 : !> \param needs the components needed by this functional are set to
83 : !> true (does not set the unneeded components to false)
84 : !> \param max_deriv ...
85 : !> \par History
86 : !> 11.2003 created [fawzi]
87 : !> \author fawzi
88 : ! **************************************************************************************************
89 2742 : SUBROUTINE lyp_lsd_info(reference, shortform, needs, max_deriv)
90 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
91 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
92 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
93 :
94 2742 : IF (PRESENT(reference)) THEN
95 42 : reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}"
96 : END IF
97 2742 : IF (PRESENT(shortform)) THEN
98 42 : shortform = "Lee-Yang-Parr correlation energy functional (LSD)"
99 : END IF
100 2742 : IF (PRESENT(needs)) THEN
101 2700 : needs%rho_spin = .TRUE.
102 2700 : needs%norm_drho_spin = .TRUE.
103 2700 : needs%norm_drho = .TRUE.
104 : END IF
105 2742 : IF (PRESENT(max_deriv)) max_deriv = 2
106 :
107 2742 : END SUBROUTINE lyp_lsd_info
108 :
109 : ! **************************************************************************************************
110 : !> \brief evaluates the lyp correlation functional for lda
111 : !> \param rho_set the density where you want to evaluate the functional
112 : !> \param deriv_set place where to store the functional derivatives (they are
113 : !> added to the derivatives)
114 : !> \param grad_deriv degree of the derivative that should be evaluated,
115 : !> if positive all the derivatives up to the given degree are evaluated,
116 : !> if negative only the given degree is calculated
117 : !> \param lyp_params input parameters (scaling)
118 : !> \par History
119 : !> 11.2003 created [fawzi]
120 : !> 01.2007 added scaling [Manuel Guidon]
121 : !> \author fawzi
122 : ! **************************************************************************************************
123 23235 : SUBROUTINE lyp_lda_eval(rho_set, deriv_set, grad_deriv, lyp_params)
124 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
125 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
126 : INTEGER, INTENT(in) :: grad_deriv
127 : TYPE(section_vals_type), POINTER :: lyp_params
128 :
129 : CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lda_eval'
130 :
131 : INTEGER :: handle, npoints
132 : INTEGER, DIMENSION(2, 3) :: bo
133 : REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho, sc
134 7745 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
135 7745 : e_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, &
136 7745 : e_rho_rho_rho, norm_drho, rho, rho_1_3
137 : TYPE(xc_derivative_type), POINTER :: deriv
138 :
139 7745 : CALL timeset(routineN, handle)
140 :
141 7745 : CALL section_vals_val_get(lyp_params, "scale_c", r_val=sc)
142 7745 : CALL cite_reference(Lee1988)
143 :
144 : CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
145 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
146 7745 : drho_cutoff=epsilon_norm_drho)
147 7745 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
148 :
149 7745 : dummy => rho
150 :
151 7745 : e_0 => dummy
152 7745 : e_rho => dummy
153 7745 : e_ndrho => dummy
154 7745 : e_rho_rho => dummy
155 7745 : e_ndrho_rho => dummy
156 7745 : e_ndrho_ndrho => dummy
157 7745 : e_rho_rho_rho => dummy
158 7745 : e_ndrho_rho_rho => dummy
159 7745 : e_ndrho_ndrho_rho => dummy
160 :
161 7745 : IF (grad_deriv >= 0) THEN
162 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
163 7745 : allocate_deriv=.TRUE.)
164 7745 : CALL xc_derivative_get(deriv, deriv_data=e_0)
165 : END IF
166 7745 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
167 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
168 7541 : allocate_deriv=.TRUE.)
169 7541 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
170 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
171 7541 : allocate_deriv=.TRUE.)
172 7541 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
173 : END IF
174 7745 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
175 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
176 390 : allocate_deriv=.TRUE.)
177 390 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
178 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
179 390 : allocate_deriv=.TRUE.)
180 390 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
181 : deriv => xc_dset_get_derivative(deriv_set, &
182 390 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
183 390 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
184 : END IF
185 7745 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
186 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
187 0 : allocate_deriv=.TRUE.)
188 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
189 : deriv => xc_dset_get_derivative(deriv_set, &
190 0 : [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
191 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
192 : deriv => xc_dset_get_derivative(deriv_set, &
193 0 : [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
194 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
195 : !FM deriv => xc_dset_get_derivative(deriv_set,&
196 : !FM [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.,&
197 : !FM
198 : !FM call xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho_ndrho)
199 : END IF
200 7745 : IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
201 0 : CPABORT("derivatives bigger than 3 not implemented")
202 : END IF
203 :
204 : !$OMP PARALLEL DEFAULT(NONE) &
205 : !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho, e_ndrho) &
206 : !$OMP SHARED(e_rho_rho, e_ndrho_rho, e_ndrho_ndrho) &
207 : !$OMP SHARED(e_rho_rho_rho, e_ndrho_rho_rho) &
208 : !$OMP SHARED(e_ndrho_ndrho_rho, grad_deriv, npoints) &
209 7745 : !$OMP SHARED(epsilon_rho, sc)
210 :
211 : CALL lyp_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
212 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
213 : e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
214 : e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
215 : e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
216 : grad_deriv=grad_deriv, &
217 : npoints=npoints, epsilon_rho=epsilon_rho, sc=sc)
218 :
219 : !$OMP END PARALLEL
220 :
221 7745 : CALL timestop(handle)
222 7745 : END SUBROUTINE lyp_lda_eval
223 :
224 : ! **************************************************************************************************
225 : !> \brief evaluates the lyp correlation functional for lda
226 : !> \param rho the density where you want to evaluate the functional
227 : !> \param rho_1_3 ...
228 : !> \param norm_drho ...
229 : !> \param e_0 ...
230 : !> \param e_rho ...
231 : !> \param e_ndrho ...
232 : !> \param e_rho_rho ...
233 : !> \param e_ndrho_rho ...
234 : !> \param e_ndrho_ndrho ...
235 : !> \param e_rho_rho_rho ...
236 : !> \param e_ndrho_rho_rho ...
237 : !> \param e_ndrho_ndrho_rho ...
238 : !> \param grad_deriv degree of the derivative that should be evaluated,
239 : !> if positive all the derivatives up to the given degree are evaluated,
240 : !> if negative only the given degree is calculated
241 : !> \param npoints ...
242 : !> \param epsilon_rho ...
243 : !> \param sc scaling-parameter for correlation
244 : !> \par History
245 : !> 11.2003 created [fawzi]
246 : !> 01.2007 added scaling [Manuel Guidon]
247 : !> \author fawzi
248 : ! **************************************************************************************************
249 7745 : SUBROUTINE lyp_lda_calc(rho, rho_1_3, norm_drho, &
250 7745 : e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
251 7745 : e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
252 : grad_deriv, npoints, epsilon_rho, sc)
253 : INTEGER, INTENT(in) :: npoints, grad_deriv
254 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_rho, e_ndrho_rho_rho, &
255 : e_rho_rho_rho, e_ndrho_ndrho, &
256 : e_ndrho_rho, e_rho_rho, e_ndrho, &
257 : e_rho, e_0
258 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
259 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sc
260 :
261 : INTEGER :: ii
262 : REAL(kind=dp) :: cf, epsilon_rho13, my_ndrho, my_rho, my_rho_1_3, t1, t102, t103, t105, &
263 : t106, t11, t112, t123, t124, t127, t128, t13, t133, t134, t14, t161, t165, t166, t173, &
264 : t176, t184, t185, t189, t192, t193, t196, t199, t2, t200, t201, t202, t203, t215, t22, &
265 : t227, t228, t231, t245, t246, t248, t26, t265, t278, t279, t3, t332, t37, t38, t39, t4, &
266 : t40, t41, t44, t45, t48, t5, t52, t6, t61, t62, t69, t7, t70, t77, t78, t80, t87, t88, &
267 : t89, t93, t94, t95, t98, t99
268 :
269 7745 : epsilon_rho13 = epsilon_rho**(1.0_dp/3.0_dp)
270 7745 : cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
271 7151 : SELECT CASE (grad_deriv)
272 : CASE (1)
273 594 : !$OMP DO
274 : DO ii = 1, npoints
275 290132085 : my_rho = rho(ii)
276 290132085 : IF (my_rho > epsilon_rho) THEN
277 277970905 : my_rho_1_3 = rho_1_3(ii)
278 277970905 : my_ndrho = norm_drho(ii)
279 277970905 : t1 = my_rho_1_3**2
280 277970905 : t2 = t1*my_rho
281 277970905 : t3 = 0.1e1_dp/t2
282 277970905 : t4 = a*t3
283 277970905 : t5 = my_rho**2
284 277970905 : t6 = t5*my_rho
285 277970905 : t7 = my_rho_1_3*t6
286 277970905 : t11 = 0.1e1_dp/my_rho_1_3
287 277970905 : t13 = EXP(-c*t11)
288 277970905 : t14 = b*t13
289 277970905 : t22 = my_ndrho**2
290 277970905 : t26 = my_rho_1_3*t22
291 : t37 = -0.72e2_dp*t7 - 0.72e2_dp*t6*d - 0.72e2_dp*t14* &
292 : t7*cf - 0.72e2_dp*t14*t6*cf*d + 0.3e1_dp*t14&
293 : & *t1*t22 + 0.10e2_dp*t14*t26*d + 0.7e1_dp*t14 &
294 277970905 : &*t26*c + 0.7e1_dp*t14*t22*c*d
295 277970905 : t38 = my_rho_1_3 + d
296 277970905 : t39 = t38**2
297 277970905 : t40 = 0.1e1_dp/t39
298 277970905 : t41 = t37*t40
299 :
300 : e_0(ii) = e_0(ii) &
301 277970905 : + (t4*t41/0.72e2_dp)*sc
302 277970905 : t44 = 0.1e1_dp/t1/t5
303 277970905 : t45 = a*t44
304 277970905 : t48 = my_rho_1_3*t5
305 277970905 : t52 = b*c
306 277970905 : t61 = t13*cf
307 277970905 : t62 = t61*d
308 277970905 : t69 = 0.1e1_dp/t1
309 277970905 : t70 = t69*t13
310 277970905 : t77 = 0.1e1_dp/my_rho
311 277970905 : t78 = t52*t77
312 277970905 : t80 = t13*t22*d
313 277970905 : t87 = c**2
314 277970905 : t88 = b*t87
315 277970905 : t89 = t77*t13
316 277970905 : t93 = my_rho_1_3*my_rho
317 277970905 : t94 = 0.1e1_dp/t93
318 277970905 : t95 = t88*t94
319 : t98 = -0.240e3_dp*t48 - 0.216e3_dp*t5*d - 0.24e2_dp*t52&
320 : & *t5*t13*cf - 0.240e3_dp*t14*t48*cf - &
321 : 0.24e2_dp*t52*t2*t62 - 0.216e3_dp*t14*t5*cf &
322 : &*d + 0.10e2_dp/0.3e1_dp*t52*t70*t22 + 0.2e1_dp* &
323 : t14*t11*t22 + 0.10e2_dp/0.3e1_dp*t78*t80 + &
324 : 0.10e2_dp/0.3e1_dp*t14*t69*t22*d + 0.7e1_dp/ &
325 : 0.3e1_dp*t88*t89*t22 + 0.7e1_dp/0.3e1_dp*t95* &
326 277970905 : t80
327 277970905 : t99 = t98*t40
328 277970905 : t102 = 0.1e1_dp/t48
329 277970905 : t103 = a*t102
330 277970905 : t105 = 0.1e1_dp/t39/t38
331 277970905 : t106 = t37*t105
332 :
333 : e_rho(ii) = e_rho(ii) &
334 : - (0.5e1_dp/0.216e3_dp*t45*t41 - t4*t99/0.72e2_dp&
335 277970905 : & + t103*t106/0.108e3_dp)*sc
336 :
337 277970905 : t112 = my_rho_1_3*my_ndrho
338 : t123 = 0.6e1_dp*t14*t1*my_ndrho + 0.20e2_dp*t14*t112 &
339 : &*d + 0.14e2_dp*t14*t112*c + 0.14e2_dp*t14* &
340 277970905 : my_ndrho*c*d
341 277970905 : t124 = t123*t40
342 :
343 : e_ndrho(ii) = e_ndrho(ii) &
344 277970905 : + (t4*t124/0.72e2_dp)*sc
345 : END IF
346 : END DO
347 : !$OMP END DO
348 : CASE default
349 7745 : !$OMP DO
350 : DO ii = 1, npoints
351 30412178 : my_rho = rho(ii)
352 30412178 : IF (my_rho > epsilon_rho) THEN
353 26486724 : my_rho_1_3 = rho_1_3(ii)
354 26486724 : my_ndrho = norm_drho(ii)
355 26486724 : t1 = my_rho_1_3**2
356 26486724 : t2 = t1*my_rho
357 26486724 : t3 = 0.1e1_dp/t2
358 26486724 : t4 = a*t3
359 26486724 : t5 = my_rho**2
360 26486724 : t6 = t5*my_rho
361 26486724 : t7 = my_rho_1_3*t6
362 26486724 : t11 = 0.1e1_dp/my_rho_1_3
363 26486724 : t13 = EXP(-c*t11)
364 26486724 : t14 = b*t13
365 26486724 : t22 = my_ndrho**2
366 26486724 : t26 = my_rho_1_3*t22
367 : t37 = -0.72e2_dp*t7 - 0.72e2_dp*t6*d - 0.72e2_dp*t14* &
368 : t7*cf - 0.72e2_dp*t14*t6*cf*d + 0.3e1_dp*t14&
369 : & *t1*t22 + 0.10e2_dp*t14*t26*d + 0.7e1_dp*t14 &
370 26486724 : &*t26*c + 0.7e1_dp*t14*t22*c*d
371 26486724 : t38 = my_rho_1_3 + d
372 26486724 : t39 = t38**2
373 26486724 : t40 = 0.1e1_dp/t39
374 26486724 : t41 = t37*t40
375 :
376 26486724 : IF (grad_deriv >= 0) THEN
377 : e_0(ii) = e_0(ii) &
378 26486724 : + (t4*t41/0.72e2_dp)*sc
379 : END IF
380 :
381 26486724 : t44 = 0.1e1_dp/t1/t5
382 26486724 : t45 = a*t44
383 26486724 : t48 = my_rho_1_3*t5
384 26486724 : t52 = b*c
385 26486724 : t61 = t13*cf
386 26486724 : t62 = t61*d
387 26486724 : t69 = 0.1e1_dp/t1
388 26486724 : t70 = t69*t13
389 26486724 : t77 = 0.1e1_dp/my_rho
390 26486724 : t78 = t52*t77
391 26486724 : t80 = t13*t22*d
392 26486724 : t87 = c**2
393 26486724 : t88 = b*t87
394 26486724 : t89 = t77*t13
395 26486724 : t93 = my_rho_1_3*my_rho
396 26486724 : t94 = 0.1e1_dp/t93
397 26486724 : t95 = t88*t94
398 : t98 = -0.240e3_dp*t48 - 0.216e3_dp*t5*d - 0.24e2_dp*t52&
399 : & *t5*t13*cf - 0.240e3_dp*t14*t48*cf - &
400 : 0.24e2_dp*t52*t2*t62 - 0.216e3_dp*t14*t5*cf &
401 : &*d + 0.10e2_dp/0.3e1_dp*t52*t70*t22 + 0.2e1_dp* &
402 : t14*t11*t22 + 0.10e2_dp/0.3e1_dp*t78*t80 + &
403 : 0.10e2_dp/0.3e1_dp*t14*t69*t22*d + 0.7e1_dp/ &
404 : 0.3e1_dp*t88*t89*t22 + 0.7e1_dp/0.3e1_dp*t95* &
405 26486724 : t80
406 26486724 : t99 = t98*t40
407 26486724 : t102 = 0.1e1_dp/t48
408 26486724 : t103 = a*t102
409 26486724 : t105 = 0.1e1_dp/t39/t38
410 26486724 : t106 = t37*t105
411 26486724 : t112 = my_rho_1_3*my_ndrho
412 : t123 = 0.6e1_dp*t14*t1*my_ndrho + 0.20e2_dp*t14*t112 &
413 : &*d + 0.14e2_dp*t14*t112*c + 0.14e2_dp*t14* &
414 26486724 : my_ndrho*c*d
415 26486724 : t124 = t123*t40
416 :
417 26486724 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
418 : e_rho(ii) = e_rho(ii) &
419 : - (0.5e1_dp/0.216e3_dp*t45*t41 - t4*t99/0.72e2_dp&
420 4320352 : & + t103*t106/0.108e3_dp)*sc
421 : e_ndrho(ii) = e_ndrho(ii) &
422 4320352 : + (t4*t124/0.72e2_dp)*sc
423 : END IF
424 :
425 26486724 : t127 = 0.1e1_dp/t1/t6
426 26486724 : t128 = a*t127
427 26486724 : t133 = 0.1e1_dp/t7
428 26486724 : t134 = a*t133
429 26486724 : t161 = t3*t13
430 26486724 : t165 = 0.1e1_dp/t5
431 26486724 : t166 = t165*t13
432 26486724 : t173 = t52*t165
433 26486724 : t176 = t88*t102
434 26486724 : t184 = b*t87*c
435 26486724 : t185 = t102*t13
436 26486724 : t189 = t184*t44
437 : t192 = -0.560e3_dp*t93 - 0.432e3_dp*my_rho*d - 0.128e3_dp&
438 : & *t52*my_rho*t13*cf - 0.8e1_dp*t88*t1*t13* &
439 : cf - 0.560e3_dp*t14*t93*cf - 0.112e3_dp*t52*t1&
440 : & *t62 - 0.8e1_dp*t88*my_rho_1_3*t62 - 0.432e3_dp* &
441 : t14*my_rho*cf*d - 0.14e2_dp/0.9e1_dp*t52* &
442 : t161*t22 - 0.11e2_dp/0.9e1_dp*t88*t166*t22 - &
443 : 0.2e1_dp/0.3e1_dp*t14*t94*t22 - 0.20e2_dp/ &
444 : 0.9e1_dp*t173*t80 - 0.2e1_dp*t176*t80 - &
445 : 0.20e2_dp/0.9e1_dp*t14*t3*t22*d + 0.7e1_dp/ &
446 : 0.9e1_dp*t184*t185*t22 + 0.7e1_dp/0.9e1_dp* &
447 26486724 : t189*t80
448 26486724 : t193 = t192*t40
449 26486724 : t196 = t98*t105
450 26486724 : t199 = 0.1e1_dp/t6
451 26486724 : t200 = a*t199
452 26486724 : t201 = t39**2
453 26486724 : t202 = 0.1e1_dp/t201
454 26486724 : t203 = t37*t202
455 26486724 : t215 = t13*my_ndrho*d
456 : t227 = 0.20e2_dp/0.3e1_dp*t52*t70*my_ndrho + 0.4e1_dp* &
457 : t14*t11*my_ndrho + 0.20e2_dp/0.3e1_dp*t78*t215&
458 : & + 0.20e2_dp/0.3e1_dp*t14*t69*my_ndrho*d + &
459 : 0.14e2_dp/0.3e1_dp*t88*t89*my_ndrho + 0.14e2_dp &
460 26486724 : &/0.3e1_dp*t95*t215
461 26486724 : t228 = t227*t40
462 26486724 : t231 = t123*t105
463 : t245 = 0.6e1_dp*t14*t1 + 0.20e2_dp*t14*my_rho_1_3*d + &
464 : 0.14e2_dp*t14*my_rho_1_3*c + 0.14e2_dp*t14*c* &
465 26486724 : d
466 26486724 : t246 = t245*t40
467 :
468 26486724 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
469 : e_rho_rho(ii) = e_rho_rho(ii) &
470 : + (0.5e1_dp/0.81e2_dp*t128*t41 - 0.5e1_dp/0.108e3_dp&
471 : & *t45*t99 + t134*t106/0.27e2_dp + t4*t193/ &
472 : 0.72e2_dp - t103*t196/0.54e2_dp + t200*t203/ &
473 4320352 : 0.108e3_dp)*sc
474 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
475 : - (0.5e1_dp/0.216e3_dp*t45*t124 - t4*t228/ &
476 4320352 : 0.72e2_dp + t103*t231/0.108e3_dp)*sc
477 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
478 4320352 : + (t4*t246/0.72e2_dp)*sc
479 : END IF
480 :
481 26486724 : t248 = t5**2
482 26486724 : t265 = 0.1e1_dp/t248
483 26486724 : t278 = t87**2
484 26486724 : t279 = b*t278
485 : t332 = -0.432e3_dp*d - 0.2240e4_dp/0.3e1_dp*my_rho_1_3 - &
486 : 0.74e2_dp/0.27e2_dp*t184*t127*t80 + 0.100e3_dp/ &
487 : 0.27e2_dp*t14*t44*t22*d + 0.7e1_dp/0.27e2_dp* &
488 : t279*t127*t13*t22 - 0.8e1_dp/0.3e1_dp*t184* &
489 : t70*cf - 0.2240e4_dp/0.3e1_dp*t14*my_rho_1_3* &
490 : cf - 0.48e2_dp*t88*t11*t13*cf - 0.944e3_dp/ &
491 : 0.3e1_dp*t52*t61 - 0.40e2_dp*t88*t69*t62 - &
492 : 0.656e3_dp/0.3e1_dp*t52*t11*t62 + 0.7e1_dp/ &
493 : 0.27e2_dp*t279*t265*t80 + 0.64e2_dp/0.27e2_dp* &
494 : t52*t44*t13*t22 - 0.8e1_dp/0.3e1_dp*t184*t77&
495 : & *t62 - 0.432e3_dp*t14*cf*d + 0.52e2_dp/ &
496 : 0.27e2_dp*t88*t199*t13*t22 - 0.20e2_dp/ &
497 : 0.9e1_dp*t184*t133*t13*t22 + 0.8e1_dp/0.9e1_dp&
498 : & *t14*t102*t22 + 0.100e3_dp/0.27e2_dp*t52*t199&
499 26486724 : & *t80 + 0.106e3_dp/0.27e2_dp*t88*t133*t80
500 :
501 26486724 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
502 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
503 : - (0.55e2_dp/0.243e3_dp*a/t1/t248*t41 - 0.5e1_dp &
504 : &/0.27e2_dp*t128*t99 + 0.40e2_dp/0.243e3_dp*a/ &
505 : my_rho_1_3/t248*t106 + 0.5e1_dp/0.72e2_dp*t45* &
506 : t193 - t134*t196/0.9e1_dp + 0.7e1_dp/0.108e3_dp* &
507 : a*t265*t203 - t4*t332*t40/0.72e2_dp + t103* &
508 : t192*t105/0.36e2_dp - t200*t98*t202/0.36e2_dp &
509 0 : &+ t128*t37/t201/t38/0.81e2_dp)*sc
510 : e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
511 : + (0.5e1_dp/0.81e2_dp*t128*t124 - 0.5e1_dp/ &
512 : 0.108e3_dp*t45*t228 + t134*t231/0.27e2_dp + t4*&
513 : & (-0.28e2_dp/0.9e1_dp*t52*t161*my_ndrho - &
514 : 0.22e2_dp/0.9e1_dp*t88*t166*my_ndrho - 0.4e1_dp &
515 : &/0.3e1_dp*t14*t94*my_ndrho - 0.40e2_dp/0.9e1_dp &
516 : &*t173*t215 - 0.4e1_dp*t176*t215 - 0.40e2_dp/ &
517 : 0.9e1_dp*t14*t3*my_ndrho*d + 0.14e2_dp/ &
518 : 0.9e1_dp*t184*t185*my_ndrho + 0.14e2_dp/0.9e1_dp&
519 : & *t189*t215)*t40/0.72e2_dp - t103*t227*t105/ &
520 0 : 0.54e2_dp + t200*t123*t202/0.108e3_dp)*sc
521 : e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
522 : - (0.5e1_dp/0.216e3_dp*t45*t246 - t4*(0.20e2_dp/ &
523 : 0.3e1_dp*t52*t70 + 0.4e1_dp*t14*t11 + 0.20e2_dp &
524 : &/0.3e1_dp*t52*t89*d + 0.20e2_dp/0.3e1_dp*t14* &
525 : t69*d + 0.14e2_dp/0.3e1_dp*t88*t89 + 0.14e2_dp/ &
526 : 0.3e1_dp*t88*t94*t13*d)*t40/0.72e2_dp + t103&
527 0 : & *t245*t105/0.108e3_dp)*sc
528 : END IF
529 : END IF
530 :
531 : !FM t1 = my_rho_1_3 ** 2
532 : !FM t2 = t1 * my_rho
533 : !FM t3 = 0.1e1_dp / t2
534 : !FM t4 = a * t3
535 : !FM t5 = my_rho ** 2
536 : !FM t6 = t5 * my_rho
537 : !FM t7 = my_rho_1_3 * t6
538 : !FM t11 = 0.1e1_dp / my_rho_1_3
539 : !FM t13 = exp(-c * t11)
540 : !FM t14 = b * t13
541 : !FM t22 = my_ndrho ** 2
542 : !FM t26 = my_rho_1_3 * t22
543 : !FM t37 = -0.72e2_dp *( t7 - t14 *&
544 : !FM & t7 * cf - t6 * d * (1.0_dp+ t14 * cf)) + t14 *(t22*(0.3e1_dp &
545 : !FM & * t1 + 0.7e1_dp * c * d)&
546 : !FM + 0.10e2_dp * t26 * d + 0.7e1_dp &
547 : !FM &* t26 * c )
548 : !FM t38 = my_rho_1_3 + d
549 : !FM t39 = t38 ** 2
550 : !FM t40 = 0.1e1_dp / t39
551 : !FM t41 = t37 * t40
552 : !FM
553 : !FM IF (grad_deriv>=0.AND.my_rho>epsilon_rho) THEN
554 : !FM e_0(ii) = e_0(ii)&
555 : !FM + t4 * t41 / 0.72e2_dp
556 : !FM END IF
557 : !FM
558 : !FM t44 = 0.1e1_dp / (t1 * t5)
559 : !FM t45 = a * t44
560 : !FM t48 = my_rho_1_3 * t5
561 : !FM t52 = b * c
562 : !FM t61 = t13 * cf
563 : !FM t62 = t61 * d
564 : !FM t69 = 0.1e1_dp / t1
565 : !FM t70 = t69 * t13
566 : !FM t77 = 0.1e1_dp / my_rho
567 : !FM t78 = t52 * t77
568 : !FM t80 = t13 * t22 * d
569 : !FM t87 = c ** 2
570 : !FM t88 = b * t87
571 : !FM t89 = t77 * t13
572 : !FM t93 = my_rho_1_3 * my_rho
573 : !FM t94 = 0.1e1_dp / t93
574 : !FM t95 = t88 * t94
575 : !FM t98 = -0.216e3_dp * t5 *d -0.240e3_dp * t48(1.0_dp+ t14 * cf) -&
576 : !FM & 0.24e2_dp * t52 * (t2 * t62 +t13*t5*cf) &
577 : !FM - 0.216e3_dp * t14 * t5 * cf &
578 : !FM &* d + t22 *(0.10e2_dp / 0.3e1_dp * t52 * t70 + 0.2e1_dp *&
579 : !FM & t14 * t11 + 0.10e2_dp / 0.3e1_dp * t14 * t69 * d + 0.7e1_dp /&
580 : !FM & 0.3e1_dp * t88 * t89 ) +(0.10e2_dp / 0.3e1_dp * t78 +&
581 : !FM & 0.7e1_dp / 0.3e1_dp * t95 ) *&
582 : !FM & t80
583 : !FM t99 = t98 * t40
584 : !FM t102 = 0.1e1_dp / t48
585 : !FM t103 = a * t102
586 : !FM t105 = 0.1e1_dp / (t39 * t38)
587 : !FM t106 = t37 * t105
588 : !FM t112 = my_rho_1_3 * my_ndrho
589 : !FM t123 = t14 *(0.6e1_dp t1 * my_ndrho + t112 * 0.20e2_dp &
590 : !FM &* d + 0.14e2_dp * c *(t112 + t14 *&
591 : !FM & my_ndrho * d))
592 : !FM t124 = t123 * t40
593 : !FM
594 : !FM IF (grad_deriv>=1.OR.grad_deriv==-1) THEN
595 : !FM e_rho(ii) = e_rho(ii) &
596 : !FM -0.5e1_dp / 0.216e3_dp * t45 * t41 + t4 * t99 / 0.72e2_dp&
597 : !FM & - t103 * t106 / 0.108e3_dp
598 : !FM e_ndrho(ii) = e_ndrho(ii)&
599 : !FM +t4 * t124 / 0.72e2_dp
600 : !FM END IF
601 : !FM
602 : !FM t127 = 0.1e1_dp / (t1 * t6)
603 : !FM t128 = a * t127
604 : !FM t133 = 0.1e1_dp / t7
605 : !FM t134 = a * t133
606 : !FM t161 = t3 * t13
607 : !FM t165 = 0.1e1_dp / t5
608 : !FM t166 = t165 * t13
609 : !FM t173 = t52 * t165
610 : !FM t176 = t88 * t102
611 : !FM t184 = b * t87 * c
612 : !FM t185 = t102 * t13
613 : !FM t189 = t184 * t44
614 : !FM t192 = -0.560e3_dp * t93 - 0.432e3_dp * my_rho * d +cf*(&
615 : !FM -t13*(0.128e3_dp&
616 : !FM & * t52 * my_rho + 0.8e1_dp * t88 * t1 )&
617 : !FM & - 0.560e3_dp * t14 * t93 ) - 0.112e3_dp * t52 * t1&
618 : !FM & * t62 - 0.8e1_dp * t88 * my_rho_1_3 * t62 - 0.432e3_dp *&
619 : !FM & t14 * my_rho * cf * d - 0.14e2_dp / 0.9e1_dp * t52 *&
620 : !FM & t161 * t22 - 0.11e2_dp / 0.9e1_dp * t88 * t166 * t22 -&
621 : !FM & 0.2e1_dp / 0.3e1_dp * t14 * t94 * t22 - 0.20e2_dp /&
622 : !FM & 0.9e1_dp * t173 * t80 - 0.2e1_dp * t176 * t80 -&
623 : !FM & 0.20e2_dp / 0.9e1_dp * t14 * t3 * t22 * d + 0.7e1_dp /&
624 : !FM & 0.9e1_dp * t184 * t185 * t22 + 0.7e1_dp / 0.9e1_dp *&
625 : !FM & t189 * t80
626 : !FM t193 = t192 * t40
627 : !FM t196 = t98 * t105
628 : !FM t199 = 0.1e1_dp / t6
629 : !FM t200 = a * t199
630 : !FM t201 = t39 ** 2
631 : !FM t202 = 0.1e1_dp / t201
632 : !FM t203 = t37 * t202
633 : !FM t215 = t13 * my_ndrho * d
634 : !FM t227 = 0.20e2_dp / 0.3e1_dp * t52 * t70 * my_ndrho + 0.4e1_dp *&
635 : !FM & t14 * t11 * my_ndrho + 0.20e2_dp / 0.3e1_dp * t78 * t215&
636 : !FM & + 0.20e2_dp / 0.3e1_dp * t14 * t69 * my_ndrho * d +&
637 : !FM & 0.14e2_dp / 0.3e1_dp * t88 * t89 * my_ndrho + 0.14e2_dp &
638 : !FM &/ 0.3e1_dp * t95 * t215
639 : !FM t228 = t227 * t40
640 : !FM t231 = t123 * t105
641 : !FM t245 = 0.6e1_dp * t14 * t1 + 0.20e2_dp * t14 * my_rho_1_3 * d +&
642 : !FM & 0.14e2_dp * t14 * my_rho_1_3 * c + 0.14e2_dp * t14 * c *&
643 : !FM & d
644 : !FM t246 = t245 * t40
645 : !FM
646 : !FM IF (grad_deriv>=2 .OR.grad_deriv==-2) THEN
647 : !FM e_rho_rho(ii) = e_rho_rho(ii)&
648 : !FM +0.5e1_dp / 0.81e2_dp * t128 * t41 - 0.5e1_dp / 0.108e3_dp&
649 : !FM & * t45 * t99 + t134 * t106 / 0.27e2_dp + t4 * t193 /&
650 : !FM & 0.72e2_dp - t103 * t196 / 0.54e2_dp + t200 * t203 /&
651 : !FM & 0.108e3_dp
652 : !FM e_ndrho_rho(ii) = e_ndrho_rho(ii)&
653 : !FM -0.5e1_dp / 0.216e3_dp * t45 * t124 + t4 * t228 /&
654 : !FM & 0.72e2_dp - t103 * t231 / 0.108e3_dp
655 : !FM e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii)&
656 : !FM +t4 * t246 / 0.72e2_dp
657 : !FM END IF
658 : !FM
659 : !FM t248 = t5 ** 2
660 : !FM t265 = 0.1e1_dp / t248
661 : !FM t278 = t87 ** 2
662 : !FM t279 = b * t278
663 : !FM t332 = -0.432e3_dp * d - 0.2240e4_dp / 0.3e1_dp * my_rho_1_3 -&
664 : !FM & 0.74e2_dp / 0.27e2_dp * t184 * t127 * t80 + 0.100e3_dp /&
665 : !FM & 0.27e2_dp * t14 * t44 * t22 * d + 0.7e1_dp / 0.27e2_dp *&
666 : !FM & t279 * t127 * t13 * t22 - 0.8e1_dp / 0.3e1_dp * t184 *&
667 : !FM & t70 * cf - 0.2240e4_dp / 0.3e1_dp * t14 * my_rho_1_3 *&
668 : !FM & cf - 0.48e2_dp * t88 * t11 * t13 * cf - 0.944e3_dp /&
669 : !FM & 0.3e1_dp * t52 * t61 - 0.40e2_dp * t88 * t69 * t62 -&
670 : !FM & 0.656e3_dp / 0.3e1_dp * t52 * t11 * t62 + 0.7e1_dp /&
671 : !FM & 0.27e2_dp * t279 * t265 * t80 + 0.64e2_dp / 0.27e2_dp *&
672 : !FM & t52 * t44 * t13 * t22 - 0.8e1_dp / 0.3e1_dp * t184 * t77&
673 : !FM & * t62 - 0.432e3_dp * t14 * cf * d + 0.52e2_dp /&
674 : !FM & 0.27e2_dp * t88 * t199 * t13 * t22 - 0.20e2_dp /&
675 : !FM & 0.9e1_dp * t184 * t133 * t13 * t22 + 0.8e1_dp / 0.9e1_dp&
676 : !FM & * t14 * t102 * t22 + 0.100e3_dp / 0.27e2_dp * t52 * t199&
677 : !FM & * t80 + 0.106e3_dp / 0.27e2_dp * t88 * t133 * t80
678 : !FM
679 : !FM IF (grad_deriv>=3 .OR. grad_deriv==-3) THEN
680 : !FM e_rho_rho_rho(ii) = e_rho_rho_rho(ii)&
681 : !FM -0.55e2_dp / 0.243e3_dp * a / t1 / t248 * t41 + 0.5e1_dp &
682 : !FM &/ 0.27e2_dp * t128 * t99 - 0.40e2_dp / 0.243e3_dp * a /&
683 : !FM & my_rho_1_3 / t248 * t106 - 0.5e1_dp / 0.72e2_dp * t45 *&
684 : !FM & t193 + t134 * t196 / 0.9e1_dp - 0.7e1_dp / 0.108e3_dp *&
685 : !FM & a * t265 * t203 + t4 * t332 * t40 / 0.72e2_dp - t103 *&
686 : !FM & t192 * t105 / 0.36e2_dp + t200 * t98 * t202 / 0.36e2_dp &
687 : !FM &- t128 * t37 / t201 / t38 / 0.81e2_dp
688 : !FM e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii)&
689 : !FM +0.5e1_dp / 0.81e2_dp * t128 * t124 - 0.5e1_dp /&
690 : !FM & 0.108e3_dp * t45 * t228 + t134 * t231 / 0.27e2_dp + t4 *&
691 : !FM & (-0.28e2_dp / 0.9e1_dp * t52 * t161 * my_ndrho -&
692 : !FM & 0.22e2_dp / 0.9e1_dp * t88 * t166 * my_ndrho - 0.4e1_dp &
693 : !FM &/ 0.3e1_dp * t14 * t94 * my_ndrho - 0.40e2_dp / 0.9e1_dp &
694 : !FM &* t173 * t215 - 0.4e1_dp * t176 * t215 - 0.40e2_dp /&
695 : !FM & 0.9e1_dp * t14 * t3 * my_ndrho * d + 0.14e2_dp /&
696 : !FM & 0.9e1_dp * t184 * t185 * my_ndrho + 0.14e2_dp / 0.9e1_dp&
697 : !FM & * t189 * t215) * t40 / 0.72e2_dp - t103 * t227 * t105 /&
698 : !FM & 0.54e2_dp + t200 * t123 * t202 / 0.108e3_dp
699 : !FM e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii)&
700 : !FM -0.5e1_dp / 0.216e3_dp * t45 * t246 + t4 * (0.20e2_dp /&
701 : !FM & 0.3e1_dp * t52 * t70 + 0.4e1_dp * t14 * t11 + 0.20e2_dp &
702 : !FM &/ 0.3e1_dp * t52 * t89 * d + 0.20e2_dp / 0.3e1_dp * t14 *&
703 : !FM & t69 * d + 0.14e2_dp / 0.3e1_dp * t88 * t89 + 0.14e2_dp /&
704 : !FM & 0.3e1_dp * t88 * t94 * t13 * d) * t40 / 0.72e2_dp - t103&
705 : !FM & * t245 * t105 / 0.108e3_dp
706 : !FM END IF
707 :
708 : END DO
709 :
710 : !$OMP END DO
711 :
712 : END SELECT
713 7745 : END SUBROUTINE lyp_lda_calc
714 :
715 : ! **************************************************************************************************
716 : !> \brief evaluates the becke 88 exchange functional for lsd
717 : !> \param rho_set ...
718 : !> \param deriv_set ...
719 : !> \param grad_deriv ...
720 : !> \param lyp_params input parameter (scaling)
721 : !> \par History
722 : !> 11.2003 created [fawzi]
723 : !> 01.2007 added scaling [Manuel Guidon]
724 : !> \author fawzi
725 : ! **************************************************************************************************
726 9042 : SUBROUTINE lyp_lsd_eval(rho_set, deriv_set, grad_deriv, lyp_params)
727 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
728 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
729 : INTEGER, INTENT(in) :: grad_deriv
730 : TYPE(section_vals_type), POINTER :: lyp_params
731 :
732 : CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lsd_eval'
733 :
734 : INTEGER :: handle, npoints
735 : INTEGER, DIMENSION(2, 3) :: bo
736 : REAL(kind=dp) :: epsilon_drho, epsilon_rho, sc
737 3014 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
738 3014 : e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, &
739 3014 : e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, &
740 3014 : norm_drhob, rhoa, rhob
741 : TYPE(xc_derivative_type), POINTER :: deriv
742 :
743 3014 : CALL timeset(routineN, handle)
744 3014 : NULLIFY (deriv)
745 :
746 3014 : CALL section_vals_val_get(lyp_params, "scale_c", r_val=sc)
747 3014 : CALL cite_reference(Lee1988)
748 :
749 : CALL xc_rho_set_get(rho_set, &
750 : rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
751 : norm_drhob=norm_drhob, norm_drho=norm_drho, &
752 : rho_cutoff=epsilon_rho, &
753 3014 : drho_cutoff=epsilon_drho, local_bounds=bo)
754 3014 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
755 :
756 3014 : dummy => rhoa
757 3014 : e_0 => dummy
758 3014 : e_ra => dummy
759 3014 : e_rb => dummy
760 3014 : e_ndra_ra => dummy
761 3014 : e_ndra_rb => dummy
762 3014 : e_ndrb_ra => dummy
763 3014 : e_ndrb_rb => dummy
764 3014 : e_ndr_ndr => dummy
765 3014 : e_ndra_ndra => dummy
766 3014 : e_ndrb_ndrb => dummy
767 3014 : e_ndr => dummy
768 3014 : e_ndra => dummy
769 3014 : e_ndrb => dummy
770 3014 : e_ra_ra => dummy
771 3014 : e_ra_rb => dummy
772 3014 : e_rb_rb => dummy
773 3014 : e_ndr_ra => dummy
774 3014 : e_ndr_rb => dummy
775 :
776 3014 : IF (grad_deriv >= 0) THEN
777 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
778 3014 : allocate_deriv=.TRUE.)
779 3014 : CALL xc_derivative_get(deriv, deriv_data=e_0)
780 : END IF
781 3014 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
782 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
783 2958 : allocate_deriv=.TRUE.)
784 2958 : CALL xc_derivative_get(deriv, deriv_data=e_ra)
785 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
786 2958 : allocate_deriv=.TRUE.)
787 2958 : CALL xc_derivative_get(deriv, deriv_data=e_rb)
788 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
789 2958 : allocate_deriv=.TRUE.)
790 2958 : CALL xc_derivative_get(deriv, deriv_data=e_ndr)
791 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
792 2958 : allocate_deriv=.TRUE.)
793 2958 : CALL xc_derivative_get(deriv, deriv_data=e_ndra)
794 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
795 2958 : allocate_deriv=.TRUE.)
796 2958 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
797 : END IF
798 3014 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
799 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
800 4 : allocate_deriv=.TRUE.)
801 4 : CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
802 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
803 4 : allocate_deriv=.TRUE.)
804 4 : CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
805 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
806 4 : allocate_deriv=.TRUE.)
807 4 : CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
808 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
809 4 : allocate_deriv=.TRUE.)
810 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndr_ra)
811 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
812 4 : allocate_deriv=.TRUE.)
813 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndr_rb)
814 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
815 4 : allocate_deriv=.TRUE.)
816 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
817 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
818 4 : allocate_deriv=.TRUE.)
819 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_rb)
820 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
821 4 : allocate_deriv=.TRUE.)
822 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
823 : deriv => xc_dset_get_derivative(deriv_set, &
824 4 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
825 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndr_ndr)
826 : deriv => xc_dset_get_derivative(deriv_set, &
827 4 : [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
828 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
829 : deriv => xc_dset_get_derivative(deriv_set, &
830 4 : [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
831 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ndrb)
832 : END IF
833 :
834 : !$OMP PARALLEL DEFAULT(NONE) &
835 : !$OMP SHARED(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob) &
836 : !$OMP SHARED(e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra) &
837 : !$OMP SHARED(e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb) &
838 : !$OMP SHARED(e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb) &
839 : !$OMP SHARED(e_ndr_ra, e_ndr_rb, grad_deriv) &
840 3014 : !$OMP SHARED(npoints, epsilon_rho, sc)
841 :
842 : CALL lyp_lsd_calc( &
843 : rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
844 : norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
845 : e_ndra_ra=e_ndra_ra, e_ndra_rb=e_ndra_rb, e_ndrb_ra&
846 : &=e_ndrb_ra, e_ndrb_rb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr, &
847 : e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr, &
848 : e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, &
849 : e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ndr_ra=e_ndr_ra, &
850 : e_ndr_rb=e_ndr_rb, &
851 : grad_deriv=grad_deriv, npoints=npoints, &
852 : epsilon_rho=epsilon_rho, sc=sc)
853 :
854 : !$OMP END PARALLEL
855 :
856 3014 : CALL timestop(handle)
857 3014 : END SUBROUTINE lyp_lsd_eval
858 :
859 : ! **************************************************************************************************
860 : !> \brief low level calculation of the becke 88 exchange functional for lsd
861 : !> \param rhoa alpha spin density
862 : !> \param rhob beta spin density
863 : !> \param norm_drho || grad rhoa + grad rhob ||
864 : !> \param norm_drhoa || grad rhoa ||
865 : !> \param norm_drhob || grad rhob ||
866 : !> \param e_0 adds to it the local value of the functional
867 : !> \param e_ra e_*: derivative of the functional wrt. to the variables
868 : !> named where the * is.
869 : !> \param e_rb ...
870 : !> \param e_ndra_ra ...
871 : !> \param e_ndra_rb ...
872 : !> \param e_ndrb_ra ...
873 : !> \param e_ndrb_rb ...
874 : !> \param e_ndr_ndr ...
875 : !> \param e_ndra_ndra ...
876 : !> \param e_ndrb_ndrb ...
877 : !> \param e_ndr ...
878 : !> \param e_ndra ...
879 : !> \param e_ndrb ...
880 : !> \param e_ra_ra ...
881 : !> \param e_ra_rb ...
882 : !> \param e_rb_rb ...
883 : !> \param e_ndr_ra ...
884 : !> \param e_ndr_rb ...
885 : !> \param grad_deriv ...
886 : !> \param npoints ...
887 : !> \param epsilon_rho ...
888 : !> \param sc scaling parameter for correlation
889 : !> \par History
890 : !> 01.2004 created [fawzi]
891 : !> \author fawzi
892 : ! **************************************************************************************************
893 3014 : SUBROUTINE lyp_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
894 : e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, &
895 : e_ndra_ndra, e_ndrb_ndrb, e_ndr, &
896 : e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb, &
897 : grad_deriv, npoints, epsilon_rho, sc)
898 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, &
899 : norm_drhob
900 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, &
901 : e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, &
902 : e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb
903 : INTEGER, INTENT(in) :: grad_deriv, npoints
904 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sc
905 :
906 : INTEGER :: ii
907 : REAL(kind=dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rho, my_rhoa, my_rhob, t1, t100, &
908 : t101, t102, t103, t104, t107, t108, t109, t112, t115, t118, t12, t120, t124, t129, t13, &
909 : t130, t132, t135, t138, t14, t140, t142, t145, t15, t153, t155, t159, t16, t164, t168, &
910 : t17, t171, t176, t18, t181, t182, t185, t189, t194, t195, t198, t2, t20, t202, t205, &
911 : t206, t21, t210, t214, t215, t218, t22, t222, t228, t23, t232, t236, t238, t24, t243, &
912 : t249, t25, t252, t253, t255, t257, t26, t260, t265, t268, t27, t270, t29, t3, t30, t304, &
913 : t31, t310, t313, t316, t319, t322, t326, t329, t332, t334, t341, t35
914 : REAL(kind=dp) :: t354, t356, t360, t363, t37, t373, t376, t381, t39, t391, t4, t40, t408, &
915 : t41, t415, t419, t422, t430, t44, t445, t449, t45, t452, t46, t467, t47, t48, t483, t487, &
916 : t49, t490, t5, t50, t505, t515, t519, t52, t520, t54, t56, t57, t6, t61, t62, t64, t66, &
917 : t7, t72, t75, t78, t8, t82, t85, t86, t87, t88, t90, t92, t94, t95, t98, t99
918 :
919 3014 : cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
920 :
921 3014 : !$OMP DO
922 :
923 : DO ii = 1, npoints
924 53011204 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
925 53011204 : my_rhob = MAX(rhob(ii), 0.0_dp)
926 53011204 : my_rho = my_rhoa + my_rhob
927 53011204 : IF (my_rho > epsilon_rho) THEN
928 51466103 : my_ndrho = norm_drho(ii)
929 51466103 : my_ndrhoa = norm_drhoa(ii)
930 51466103 : my_ndrhob = norm_drhob(ii)
931 :
932 51466103 : t1 = my_rho**(0.1e1_dp/0.3e1_dp)
933 51466103 : t2 = 0.1e1_dp/t1
934 51466103 : t3 = d*t2
935 51466103 : t4 = 0.1e1_dp + t3
936 51466103 : t5 = 0.1e1_dp/t4
937 51466103 : t6 = a*t5
938 51466103 : t7 = my_rhoa*my_rhob
939 51466103 : t8 = 0.1e1_dp/my_rho
940 51466103 : t12 = a*b
941 51466103 : t13 = c*t2
942 51466103 : t14 = EXP(-t13)
943 51466103 : t15 = t12*t14
944 51466103 : t16 = my_rho**2
945 51466103 : t17 = t16*my_rho
946 51466103 : t18 = t1**2
947 51466103 : t20 = 0.1e1_dp/t18/t17
948 51466103 : t21 = t5*t20
949 51466103 : t22 = 2**(0.1e1_dp/0.3e1_dp)
950 51466103 : t23 = t22**2
951 51466103 : t24 = t23*cf
952 51466103 : t25 = my_rhoa**2
953 51466103 : t26 = my_rhoa**(0.1e1_dp/0.3e1_dp)
954 51466103 : t27 = t26**2
955 51466103 : t29 = my_rhob**2
956 51466103 : t30 = my_rhob**(0.1e1_dp/0.3e1_dp)
957 51466103 : t31 = t30**2
958 51466103 : t35 = 0.8e1_dp*t24*(t27*t25 + t31*t29)
959 51466103 : t37 = t3*t5
960 : t39 = 0.47e2_dp/0.18e2_dp - 0.7e1_dp/0.18e2_dp*t13 - &
961 51466103 : 0.7e1_dp/0.18e2_dp*t37
962 51466103 : t40 = my_ndrho**2
963 51466103 : t41 = t39*t40
964 51466103 : t44 = 0.5e1_dp/0.2e1_dp - t13/0.18e2_dp - t37/0.18e2_dp
965 51466103 : t45 = my_ndrhoa**2
966 51466103 : t46 = my_ndrhob**2
967 51466103 : t47 = t45 + t46
968 51466103 : t48 = t44*t47
969 51466103 : t49 = t13 + t37 - 0.11e2_dp
970 51466103 : t50 = my_rhoa*t8
971 51466103 : t52 = my_rhob*t8
972 51466103 : t54 = t50*t45 + t52*t46
973 51466103 : t56 = t49*t54/0.9e1_dp
974 51466103 : t57 = t35 + t41 - t48 - t56
975 51466103 : t61 = 0.2e1_dp/0.3e1_dp*t16
976 51466103 : t62 = t61 - t25
977 51466103 : t64 = t61 - t29
978 : t66 = t7*t57 - 0.2e1_dp/0.3e1_dp*t16*t40 + t62*t46 + &
979 51466103 : t64*t45
980 :
981 51466103 : IF (grad_deriv >= 0 .AND. my_rho > epsilon_rho) THEN
982 : e_0(ii) = e_0(ii) &
983 51466103 : - (0.4e1_dp*t6*t7*t8 + t15*t21*t66)*sc
984 : END IF
985 : !--------
986 :
987 51466103 : t72 = t27*my_rhoa
988 51466103 : t75 = t49*t8
989 51466103 : t78 = 0.64e2_dp/0.3e1_dp*t24*t72 - t75*t45/0.9e1_dp
990 51466103 : t82 = my_rhob*t57 + t7*t78 - 0.2e1_dp*my_rhoa*t46
991 51466103 : t85 = t4**2
992 51466103 : t86 = 0.1e1_dp/t85
993 51466103 : t87 = a*t86
994 51466103 : t88 = t87*my_rhoa
995 51466103 : t90 = 0.1e1_dp/t1/t16
996 51466103 : t92 = my_rhob*t90*d
997 51466103 : t94 = 0.4e1_dp/0.3e1_dp*t88*t92
998 51466103 : t95 = 0.1e1_dp/t16
999 51466103 : t98 = 0.4e1_dp*t6*t7*t95
1000 51466103 : t99 = t12*c
1001 51466103 : t100 = t16**2
1002 51466103 : t101 = t100*my_rho
1003 51466103 : t102 = 0.1e1_dp/t101
1004 51466103 : t103 = t102*t14
1005 51466103 : t104 = t5*t66
1006 51466103 : t107 = t99*t103*t104/0.3e1_dp
1007 51466103 : t108 = t86*t102
1008 51466103 : t109 = t66*d
1009 51466103 : t112 = t15*t108*t109/0.3e1_dp
1010 51466103 : t115 = t5/t18/t100
1011 51466103 : t118 = 0.11e2_dp/0.3e1_dp*t15*t115*t66
1012 51466103 : t120 = 0.1e1_dp/t1/my_rho
1013 51466103 : t124 = d**2
1014 51466103 : t129 = c*t120 + d*t120*t5 - t124/t18/my_rho*t86
1015 51466103 : t130 = 0.7e1_dp/0.54e2_dp*t129
1016 51466103 : t132 = t129/0.54e2_dp
1017 51466103 : t135 = -t129/0.3e1_dp
1018 51466103 : t138 = my_rhoa*t95
1019 51466103 : t140 = my_rhob*t95
1020 51466103 : t142 = -t138*t45 - t140*t46
1021 : t145 = t130*t40 - t132*t47 - t135*t54/0.9e1_dp - t49* &
1022 51466103 : t142/0.9e1_dp
1023 : t153 = t7*t145 - 0.4e1_dp/0.3e1_dp*my_rho*t40 + &
1024 : 0.4e1_dp/0.3e1_dp*my_rho*t46 + 0.4e1_dp/0.3e1_dp&
1025 51466103 : & *my_rho*t45
1026 51466103 : t155 = t15*t21*t153
1027 :
1028 51466103 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1029 : e_ra(ii) = e_ra(ii) &
1030 : - (0.4e1_dp*t6*t52 + t15*t21*t82 + t94 - t98 + &
1031 47510935 : t107 + t112 - t118 + t155)*sc
1032 : END IF
1033 :
1034 51466103 : t159 = t31*my_rhob
1035 51466103 : t164 = 0.64e2_dp/0.3e1_dp*t24*t159 - t75*t46/0.9e1_dp
1036 51466103 : t168 = my_rhoa*t57 + t7*t164 - 0.2e1_dp*my_rhob*t45
1037 :
1038 51466103 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1039 : e_rb(ii) = e_rb(ii) &
1040 : - (0.4e1_dp*t6*t50 + t15*t21*t168 + t94 - t98 + &
1041 47510935 : t107 + t112 - t118 + t155)*sc
1042 : END IF
1043 :
1044 51466103 : t171 = t39*my_ndrho
1045 51466103 : t176 = 0.2e1_dp*t7*t171 - 0.4e1_dp/0.3e1_dp*t16*my_ndrho
1046 :
1047 51466103 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1048 : e_ndr(ii) = e_ndr(ii) &
1049 47510935 : - (t15*t21*t176)*sc
1050 : END IF
1051 :
1052 51466103 : t181 = t49*my_rhoa
1053 51466103 : t182 = t8*my_ndrhoa
1054 51466103 : t185 = -0.2e1_dp*t44*my_ndrhoa - 0.2e1_dp/0.9e1_dp*t181*t182
1055 51466103 : t189 = t7*t185 + 0.2e1_dp*t64*my_ndrhoa
1056 :
1057 51466103 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1058 : e_ndra(ii) = e_ndra(ii) &
1059 47510935 : - (t15*t21*t189)*sc
1060 : END IF
1061 :
1062 51466103 : t194 = t49*my_rhob
1063 51466103 : t195 = t8*my_ndrhob
1064 51466103 : t198 = -0.2e1_dp*t44*my_ndrhob - 0.2e1_dp/0.9e1_dp*t194*t195
1065 51466103 : t202 = t7*t198 + 0.2e1_dp*t62*my_ndrhob
1066 :
1067 51466103 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
1068 : e_ndrb(ii) = e_ndrb(ii) &
1069 47510935 : - (t15*t21*t202)*sc
1070 : END IF
1071 : !-------
1072 :
1073 51466103 : t205 = t100*t16
1074 51466103 : t206 = 0.1e1_dp/t205
1075 51466103 : t210 = 0.26e2_dp/0.9e1_dp*t99*t206*t14*t104
1076 51466103 : t214 = 0.2e1_dp/0.3e1_dp*t99*t103*t5*t153
1077 51466103 : t215 = c**2
1078 51466103 : t218 = 0.1e1_dp/t1/t205
1079 51466103 : t222 = t12*t215*t218*t14*t104/0.9e1_dp
1080 51466103 : t228 = 0.2e1_dp/0.9e1_dp*t12*c*t218*t14*t86*t109
1081 51466103 : t232 = 0.26e2_dp/0.9e1_dp*t15*t86*t206*t109
1082 51466103 : t236 = 0.2e1_dp/0.3e1_dp*t15*t108*t153*d
1083 51466103 : t238 = 0.1e1_dp/t85/t4
1084 51466103 : t243 = 0.2e1_dp/0.9e1_dp*t15*t238*t218*t66*t124
1085 51466103 : t249 = 0.154e3_dp/0.9e1_dp*t15*t5/t18/t101*t66
1086 51466103 : t252 = 0.22e2_dp/0.3e1_dp*t15*t115*t153
1087 51466103 : t253 = t6*t140
1088 51466103 : t255 = t87*t92
1089 51466103 : t257 = c*t90
1090 51466103 : t260 = d*t90*t5
1091 51466103 : t265 = t124/t18/t16*t86
1092 51466103 : t268 = 0.1e1_dp/t17
1093 51466103 : t270 = t124*d*t268*t238
1094 : t304 = t15*t21*(t7*((-0.14e2_dp/0.81e2_dp*t257 - &
1095 : 0.14e2_dp/0.81e2_dp*t260 + 0.7e1_dp/0.27e2_dp* &
1096 : t265 - 0.7e1_dp/0.81e2_dp*t270)*t40 - (-0.2e1_dp/ &
1097 : 0.81e2_dp*t257 - 0.2e1_dp/0.81e2_dp*t260 + t265/ &
1098 : 0.27e2_dp - t270/0.81e2_dp)*t47 - (0.4e1_dp/ &
1099 : 0.9e1_dp*t257 + 0.4e1_dp/0.9e1_dp*t260 - 0.2e1_dp &
1100 : &/0.3e1_dp*t265 + 0.2e1_dp/0.9e1_dp*t270)*t54/ &
1101 : 0.9e1_dp - 0.2e1_dp/0.9e1_dp*t135*t142 - t49*&
1102 : & (0.2e1_dp*my_rhoa*t268*t45 + 0.2e1_dp*my_rhob* &
1103 : t268*t46)/0.9e1_dp) - 0.4e1_dp/0.3e1_dp*t40 + &
1104 51466103 : 0.4e1_dp/0.3e1_dp*t46 + 0.4e1_dp/0.3e1_dp*t45)
1105 51466103 : t310 = 0.40e2_dp/0.9e1_dp*t88*my_rhob/t1/t17*d
1106 51466103 : t313 = my_rhob*t20
1107 51466103 : t316 = 0.8e1_dp/0.9e1_dp*a*t238*my_rhoa*t313*t124
1108 51466103 : t319 = 0.8e1_dp*t6*t7*t268
1109 51466103 : t322 = t99*t103*t5*t82
1110 51466103 : t326 = t15*t108*t82*d
1111 51466103 : t329 = t15*t115*t82
1112 51466103 : t332 = t135*t8
1113 51466103 : t334 = t49*t95
1114 : t341 = t15*t21*(my_rhob*t145 + t7*(-t332*t45/ &
1115 51466103 : 0.9e1_dp + t334*t45/0.9e1_dp))
1116 :
1117 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1118 : e_ra_ra(ii) = e_ra_ra(ii) &
1119 : + (t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 + &
1120 : t252 + 0.8e1_dp*t253 - 0.8e1_dp/0.3e1_dp*t255 - &
1121 : t304 + t310 - t316 - t319 - 0.2e1_dp/0.3e1_dp*t322 - &
1122 : 0.2e1_dp/0.3e1_dp*t326 + 0.22e2_dp/0.3e1_dp*t329&
1123 : & - 0.2e1_dp*t341 - t15*t21*(0.2e1_dp*my_rhob* &
1124 : t78 + 0.320e3_dp/0.9e1_dp*t72*my_rhob*t24 - &
1125 79625 : 0.2e1_dp*t46))*sc
1126 : END IF
1127 :
1128 51466103 : t354 = t99*t103*t5*t168
1129 51466103 : t356 = t6*t138
1130 51466103 : t360 = t87*my_rhoa*t90*d
1131 51466103 : t363 = t15*t115*t168
1132 : t373 = t15*t21*(my_rhoa*t145 + t7*(-t332*t46/ &
1133 51466103 : 0.9e1_dp + t334*t46/0.9e1_dp))
1134 51466103 : t376 = t15*t108*t168*d
1135 :
1136 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1137 : e_rb_rb(ii) = e_rb_rb(ii) &
1138 : + (t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 + &
1139 : t252 - t304 + t310 - t316 - t319 + 0.8e1_dp*t356 - &
1140 : 0.8e1_dp/0.3e1_dp*t360 + 0.22e2_dp/0.3e1_dp* &
1141 : t363 - 0.2e1_dp*t373 - 0.2e1_dp/0.3e1_dp*t354 - &
1142 : 0.2e1_dp/0.3e1_dp*t376 - t15*t21*(0.2e1_dp* &
1143 : my_rhoa*t164 + 0.320e3_dp/0.9e1_dp*my_rhoa* &
1144 79625 : t159*t24 - 0.2e1_dp*t45))*sc
1145 : END IF
1146 :
1147 : t381 = -t354/0.3e1_dp + 0.4e1_dp*t356 - 0.4e1_dp/0.3e1_dp&
1148 : & *t360 + 0.11e2_dp/0.3e1_dp*t363 - t373 - t341 - &
1149 : t376/0.3e1_dp + t310 - 0.4e1_dp*t6*t8 + 0.4e1_dp* &
1150 51466103 : t253 + t210 - t214 - t222
1151 : t391 = -t228 + t232 - t236 - t243 - t249 + t252 - 0.4e1_dp/ &
1152 : 0.3e1_dp*t255 - t304 - t316 - t319 - t322/0.3e1_dp - &
1153 : t326/0.3e1_dp + 0.11e2_dp/0.3e1_dp*t329 - t15* &
1154 : t21*(t35 + t41 - t48 - t56 + my_rhob*t164 + my_rhoa &
1155 51466103 : &*t78)
1156 :
1157 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1158 : e_ra_rb(ii) = e_ra_rb(ii) &
1159 79625 : + (t381 + t391)*sc
1160 : END IF
1161 :
1162 51466103 : t408 = t12*t14*t5
1163 51466103 : t415 = t99*t103*t5*t176/0.3e1_dp
1164 51466103 : t419 = t15*t108*t176*d/0.3e1_dp
1165 51466103 : t422 = 0.11e2_dp/0.3e1_dp*t15*t115*t176
1166 : t430 = t15*t21*(0.2e1_dp*t7*t130*my_ndrho - 0.8e1_dp &
1167 51466103 : &/0.3e1_dp*my_rho*my_ndrho)
1168 :
1169 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1170 : e_ndr_ra(ii) = e_ndr_ra(ii) &
1171 : - (0.2e1_dp*t408*t313*t171 + t415 + t419 - t422 + &
1172 79625 : t430)*sc
1173 : e_ndr_rb(ii) = e_ndr_rb(ii) &
1174 : - (0.2e1_dp*t408*t20*my_rhoa*t171 + t415 + t419 - &
1175 79625 : t422 + t430)*sc
1176 : END IF
1177 :
1178 51466103 : t445 = t99*t103*t5*t189/0.3e1_dp
1179 51466103 : t449 = t15*t108*t189*d/0.3e1_dp
1180 51466103 : t452 = 0.11e2_dp/0.3e1_dp*t15*t115*t189
1181 : t467 = t15*t21*(t7*(-0.2e1_dp*t132*my_ndrhoa - &
1182 : 0.2e1_dp/0.9e1_dp*t135*my_rhoa*t182 + 0.2e1_dp/ &
1183 : 0.9e1_dp*t181*t95*my_ndrhoa) + 0.8e1_dp/0.3e1_dp&
1184 51466103 : & *my_rho*my_ndrhoa)
1185 :
1186 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1187 : e_ndra_ra(ii) = e_ndra_ra(ii) &
1188 : - (t15*t21*(my_rhob*t185 - 0.2e1_dp/0.9e1_dp*t7&
1189 79625 : & *t75*my_ndrhoa) + t445 + t449 - t452 + t467)*sc
1190 : e_ndra_rb(ii) = e_ndra_rb(ii) &
1191 : - (t15*t21*(my_rhoa*t185 - 0.4e1_dp*my_rhob* &
1192 79625 : my_ndrhoa) + t445 + t449 - t452 + t467)*sc
1193 : END IF
1194 :
1195 51466103 : t483 = t99*t103*t5*t202/0.3e1_dp
1196 51466103 : t487 = t15*t108*t202*d/0.3e1_dp
1197 51466103 : t490 = 0.11e2_dp/0.3e1_dp*t15*t115*t202
1198 : t505 = t15*t21*(t7*(-0.2e1_dp*t132*my_ndrhob - &
1199 : 0.2e1_dp/0.9e1_dp*t135*my_rhob*t195 + 0.2e1_dp/ &
1200 : 0.9e1_dp*t194*t95*my_ndrhob) + 0.8e1_dp/0.3e1_dp&
1201 51466103 : & *my_rho*my_ndrhob)
1202 :
1203 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1204 : e_ndrb_ra(ii) = e_ndrb_ra(ii) &
1205 : - (t15*t21*(my_rhob*t198 - 0.4e1_dp*my_rhoa* &
1206 79625 : my_ndrhob) + t483 + t487 - t490 + t505)*sc
1207 : e_ndrb_rb(ii) = e_ndrb_rb(ii) &
1208 : - (t15*t21*(my_rhoa*t198 - 0.2e1_dp/0.9e1_dp*t7&
1209 79625 : & *t75*my_ndrhob) + t483 + t487 - t490 + t505)*sc
1210 : END IF
1211 :
1212 51466103 : t515 = 0.4e1_dp/0.3e1_dp*t16
1213 :
1214 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1215 : e_ndr_ndr(ii) = e_ndr_ndr(ii) &
1216 79625 : - (t15*t21*(0.2e1_dp*t7*t39 - t515))*sc
1217 : END IF
1218 :
1219 51466103 : t519 = t13/0.9e1_dp
1220 51466103 : t520 = t37/0.9e1_dp
1221 :
1222 51466103 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
1223 : e_ndra_ndra(ii) = e_ndra_ndra(ii) &
1224 : - (t15*t21*(t7*(-0.5e1_dp + t519 + t520 - 0.2e1_dp &
1225 79625 : &/0.9e1_dp*t181*t8) + t515 - 0.2e1_dp*t29))*sc
1226 :
1227 : e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) &
1228 : - (t15*t21*(t7*(-0.5e1_dp + t519 + t520 - 0.2e1_dp &
1229 79625 : &/0.9e1_dp*t194*t8) + t515 - 0.2e1_dp*t25))*sc
1230 : END IF
1231 : END IF
1232 : END DO
1233 : !$OMP END DO
1234 :
1235 3014 : END SUBROUTINE lyp_lsd_calc
1236 :
1237 : END MODULE xc_lyp
|
