Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 8620 : 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 8620 : IF (PRESENT(reference)) THEN
64 105 : reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}"
65 : END IF
66 8620 : IF (PRESENT(shortform)) THEN
67 105 : shortform = "Lee-Yang-Parr correlation energy functional (LDA)"
68 : END IF
69 8620 : IF (PRESENT(needs)) THEN
70 8515 : needs%rho = .TRUE.
71 8515 : needs%rho_1_3 = .TRUE.
72 8515 : needs%norm_drho = .TRUE.
73 : END IF
74 8620 : IF (PRESENT(max_deriv)) max_deriv = 3
75 :
76 8620 : 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 2826 : 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 2826 : IF (PRESENT(reference)) THEN
95 40 : reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}"
96 : END IF
97 2826 : IF (PRESENT(shortform)) THEN
98 40 : shortform = "Lee-Yang-Parr correlation energy functional (LSD)"
99 : END IF
100 2826 : IF (PRESENT(needs)) THEN
101 2786 : needs%rho_spin = .TRUE.
102 2786 : needs%norm_drho_spin = .TRUE.
103 2786 : needs%norm_drho = .TRUE.
104 : END IF
105 2826 : IF (PRESENT(max_deriv)) max_deriv = 2
106 :
107 2826 : 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 maximum degree of the derivative that should be evaluated
115 : !> \param lyp_params input parameters (scaling)
116 : !> \par History
117 : !> 11.2003 created [fawzi]
118 : !> 01.2007 added scaling [Manuel Guidon]
119 : !> \author fawzi
120 : ! **************************************************************************************************
121 24831 : SUBROUTINE lyp_lda_eval(rho_set, deriv_set, grad_deriv, lyp_params)
122 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
123 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
124 : INTEGER, INTENT(in) :: grad_deriv
125 : TYPE(section_vals_type), POINTER :: lyp_params
126 :
127 : CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lda_eval'
128 :
129 : INTEGER :: handle, npoints
130 : INTEGER, DIMENSION(2, 3) :: bo
131 : REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho, sc
132 8277 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
133 8277 : e_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, &
134 8277 : e_rho_rho_rho, norm_drho, rho, rho_1_3
135 : TYPE(xc_derivative_type), POINTER :: deriv
136 :
137 8277 : CALL timeset(routineN, handle)
138 :
139 8277 : CALL section_vals_val_get(lyp_params, "scale_c", r_val=sc)
140 8277 : CALL cite_reference(Lee1988)
141 :
142 : CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
143 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
144 8277 : drho_cutoff=epsilon_norm_drho)
145 8277 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
146 :
147 8277 : dummy => rho
148 :
149 8277 : e_0 => dummy
150 8277 : e_rho => dummy
151 8277 : e_ndrho => dummy
152 8277 : e_rho_rho => dummy
153 8277 : e_ndrho_rho => dummy
154 8277 : e_ndrho_ndrho => dummy
155 8277 : e_rho_rho_rho => dummy
156 8277 : e_ndrho_rho_rho => dummy
157 8277 : e_ndrho_ndrho_rho => dummy
158 :
159 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
160 8277 : allocate_deriv=.TRUE.)
161 8277 : CALL xc_derivative_get(deriv, deriv_data=e_0)
162 8277 : IF (grad_deriv >= 1) THEN
163 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
164 8059 : allocate_deriv=.TRUE.)
165 8059 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
166 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
167 8059 : allocate_deriv=.TRUE.)
168 8059 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
169 : END IF
170 8277 : IF (grad_deriv >= 2) THEN
171 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
172 276 : allocate_deriv=.TRUE.)
173 276 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
174 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
175 276 : allocate_deriv=.TRUE.)
176 276 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
177 : deriv => xc_dset_get_derivative(deriv_set, &
178 276 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
179 276 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
180 : END IF
181 8277 : IF (grad_deriv >= 3) THEN
182 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
183 0 : allocate_deriv=.TRUE.)
184 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
185 : deriv => xc_dset_get_derivative(deriv_set, &
186 0 : [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
187 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
188 : deriv => xc_dset_get_derivative(deriv_set, &
189 0 : [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
190 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
191 : !FM deriv => xc_dset_get_derivative(deriv_set,&
192 : !FM [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.,&
193 : !FM
194 : !FM call xc_derivative_get(deriv,deriv_data=e_ndrho_ndrho_ndrho)
195 : END IF
196 8277 : IF (grad_deriv > 3) THEN
197 0 : CPABORT("derivatives bigger than 3 not implemented")
198 : END IF
199 :
200 : !$OMP PARALLEL DEFAULT(NONE) &
201 : !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho, e_ndrho) &
202 : !$OMP SHARED(e_rho_rho, e_ndrho_rho, e_ndrho_ndrho) &
203 : !$OMP SHARED(e_rho_rho_rho, e_ndrho_rho_rho) &
204 : !$OMP SHARED(e_ndrho_ndrho_rho, grad_deriv, npoints) &
205 8277 : !$OMP SHARED(epsilon_rho, sc)
206 :
207 : CALL lyp_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
208 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
209 : e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
210 : e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
211 : e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
212 : grad_deriv=grad_deriv, &
213 : npoints=npoints, epsilon_rho=epsilon_rho, sc=sc)
214 :
215 : !$OMP END PARALLEL
216 :
217 8277 : CALL timestop(handle)
218 8277 : END SUBROUTINE lyp_lda_eval
219 :
220 : ! **************************************************************************************************
221 : !> \brief evaluates the lyp correlation functional for lda
222 : !> \param rho the density where you want to evaluate the functional
223 : !> \param rho_1_3 ...
224 : !> \param norm_drho ...
225 : !> \param e_0 ...
226 : !> \param e_rho ...
227 : !> \param e_ndrho ...
228 : !> \param e_rho_rho ...
229 : !> \param e_ndrho_rho ...
230 : !> \param e_ndrho_ndrho ...
231 : !> \param e_rho_rho_rho ...
232 : !> \param e_ndrho_rho_rho ...
233 : !> \param e_ndrho_ndrho_rho ...
234 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
235 : !> \param npoints ...
236 : !> \param epsilon_rho ...
237 : !> \param sc scaling-parameter for correlation
238 : !> \par History
239 : !> 11.2003 created [fawzi]
240 : !> 01.2007 added scaling [Manuel Guidon]
241 : !> \author fawzi
242 : ! **************************************************************************************************
243 8277 : SUBROUTINE lyp_lda_calc(rho, rho_1_3, norm_drho, &
244 8277 : e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
245 8277 : e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
246 : grad_deriv, npoints, epsilon_rho, sc)
247 : INTEGER, INTENT(in) :: npoints, grad_deriv
248 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_rho, e_ndrho_rho_rho, &
249 : e_rho_rho_rho, e_ndrho_ndrho, &
250 : e_ndrho_rho, e_rho_rho, e_ndrho, &
251 : e_rho, e_0
252 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
253 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sc
254 :
255 : INTEGER :: ii
256 : REAL(kind=dp) :: cf, epsilon_rho13, my_ndrho, my_rho, my_rho_1_3, t1, t102, t103, t105, &
257 : t106, t11, t112, t123, t124, t127, t128, t13, t133, t134, t14, t161, t165, t166, t173, &
258 : t176, t184, t185, t189, t192, t193, t196, t199, t2, t200, t201, t202, t203, t215, t22, &
259 : t227, t228, t231, t245, t246, t248, t26, t265, t278, t279, t3, t332, t37, t38, t39, t4, &
260 : t40, t41, t44, t45, t48, t5, t52, t6, t61, t62, t69, t7, t70, t77, t78, t80, t87, t88, &
261 : t89, t93, t94, t95, t98, t99
262 :
263 8277 : epsilon_rho13 = epsilon_rho**(1.0_dp/3.0_dp)
264 8277 : cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
265 7783 : SELECT CASE (grad_deriv)
266 : CASE (1)
267 494 : !$OMP DO
268 : DO ii = 1, npoints
269 390747348 : my_rho = rho(ii)
270 390747348 : IF (my_rho > epsilon_rho) THEN
271 375880352 : my_rho_1_3 = rho_1_3(ii)
272 375880352 : my_ndrho = norm_drho(ii)
273 375880352 : t1 = my_rho_1_3**2
274 375880352 : t2 = t1*my_rho
275 375880352 : t3 = 0.1e1_dp/t2
276 375880352 : t4 = a*t3
277 375880352 : t5 = my_rho**2
278 375880352 : t6 = t5*my_rho
279 375880352 : t7 = my_rho_1_3*t6
280 375880352 : t11 = 0.1e1_dp/my_rho_1_3
281 375880352 : t13 = EXP(-c*t11)
282 375880352 : t14 = b*t13
283 375880352 : t22 = my_ndrho**2
284 375880352 : t26 = my_rho_1_3*t22
285 : t37 = -0.72e2_dp*t7 - 0.72e2_dp*t6*d - 0.72e2_dp*t14* &
286 : t7*cf - 0.72e2_dp*t14*t6*cf*d + 0.3e1_dp*t14&
287 : & *t1*t22 + 0.10e2_dp*t14*t26*d + 0.7e1_dp*t14 &
288 375880352 : &*t26*c + 0.7e1_dp*t14*t22*c*d
289 375880352 : t38 = my_rho_1_3 + d
290 375880352 : t39 = t38**2
291 375880352 : t40 = 0.1e1_dp/t39
292 375880352 : t41 = t37*t40
293 :
294 : e_0(ii) = e_0(ii) &
295 375880352 : + (t4*t41/0.72e2_dp)*sc
296 375880352 : t44 = 0.1e1_dp/t1/t5
297 375880352 : t45 = a*t44
298 375880352 : t48 = my_rho_1_3*t5
299 375880352 : t52 = b*c
300 375880352 : t61 = t13*cf
301 375880352 : t62 = t61*d
302 375880352 : t69 = 0.1e1_dp/t1
303 375880352 : t70 = t69*t13
304 375880352 : t77 = 0.1e1_dp/my_rho
305 375880352 : t78 = t52*t77
306 375880352 : t80 = t13*t22*d
307 375880352 : t87 = c**2
308 375880352 : t88 = b*t87
309 375880352 : t89 = t77*t13
310 375880352 : t93 = my_rho_1_3*my_rho
311 375880352 : t94 = 0.1e1_dp/t93
312 375880352 : t95 = t88*t94
313 : t98 = -0.240e3_dp*t48 - 0.216e3_dp*t5*d - 0.24e2_dp*t52&
314 : & *t5*t13*cf - 0.240e3_dp*t14*t48*cf - &
315 : 0.24e2_dp*t52*t2*t62 - 0.216e3_dp*t14*t5*cf &
316 : &*d + 0.10e2_dp/0.3e1_dp*t52*t70*t22 + 0.2e1_dp* &
317 : t14*t11*t22 + 0.10e2_dp/0.3e1_dp*t78*t80 + &
318 : 0.10e2_dp/0.3e1_dp*t14*t69*t22*d + 0.7e1_dp/ &
319 : 0.3e1_dp*t88*t89*t22 + 0.7e1_dp/0.3e1_dp*t95* &
320 375880352 : t80
321 375880352 : t99 = t98*t40
322 375880352 : t102 = 0.1e1_dp/t48
323 375880352 : t103 = a*t102
324 375880352 : t105 = 0.1e1_dp/t39/t38
325 375880352 : t106 = t37*t105
326 :
327 : e_rho(ii) = e_rho(ii) &
328 : - (0.5e1_dp/0.216e3_dp*t45*t41 - t4*t99/0.72e2_dp&
329 375880352 : & + t103*t106/0.108e3_dp)*sc
330 :
331 375880352 : t112 = my_rho_1_3*my_ndrho
332 : t123 = 0.6e1_dp*t14*t1*my_ndrho + 0.20e2_dp*t14*t112 &
333 : &*d + 0.14e2_dp*t14*t112*c + 0.14e2_dp*t14* &
334 375880352 : my_ndrho*c*d
335 375880352 : t124 = t123*t40
336 :
337 : e_ndrho(ii) = e_ndrho(ii) &
338 375880352 : + (t4*t124/0.72e2_dp)*sc
339 : END IF
340 : END DO
341 : !$OMP END DO
342 : CASE default
343 8277 : !$OMP DO
344 : DO ii = 1, npoints
345 30840674 : my_rho = rho(ii)
346 30840674 : IF (my_rho > epsilon_rho) THEN
347 26338552 : my_rho_1_3 = rho_1_3(ii)
348 26338552 : my_ndrho = norm_drho(ii)
349 26338552 : t1 = my_rho_1_3**2
350 26338552 : t2 = t1*my_rho
351 26338552 : t3 = 0.1e1_dp/t2
352 26338552 : t4 = a*t3
353 26338552 : t5 = my_rho**2
354 26338552 : t6 = t5*my_rho
355 26338552 : t7 = my_rho_1_3*t6
356 26338552 : t11 = 0.1e1_dp/my_rho_1_3
357 26338552 : t13 = EXP(-c*t11)
358 26338552 : t14 = b*t13
359 26338552 : t22 = my_ndrho**2
360 26338552 : t26 = my_rho_1_3*t22
361 : t37 = -0.72e2_dp*t7 - 0.72e2_dp*t6*d - 0.72e2_dp*t14* &
362 : t7*cf - 0.72e2_dp*t14*t6*cf*d + 0.3e1_dp*t14&
363 : & *t1*t22 + 0.10e2_dp*t14*t26*d + 0.7e1_dp*t14 &
364 26338552 : &*t26*c + 0.7e1_dp*t14*t22*c*d
365 26338552 : t38 = my_rho_1_3 + d
366 26338552 : t39 = t38**2
367 26338552 : t40 = 0.1e1_dp/t39
368 26338552 : t41 = t37*t40
369 :
370 : e_0(ii) = e_0(ii) &
371 26338552 : + (t4*t41/0.72e2_dp)*sc
372 :
373 26338552 : IF (grad_deriv > 0) THEN
374 2297800 : t44 = 0.1e1_dp/t1/t5
375 2297800 : t45 = a*t44
376 2297800 : t48 = my_rho_1_3*t5
377 2297800 : t52 = b*c
378 2297800 : t61 = t13*cf
379 2297800 : t62 = t61*d
380 2297800 : t69 = 0.1e1_dp/t1
381 2297800 : t70 = t69*t13
382 2297800 : t77 = 0.1e1_dp/my_rho
383 2297800 : t78 = t52*t77
384 2297800 : t80 = t13*t22*d
385 2297800 : t87 = c**2
386 2297800 : t88 = b*t87
387 2297800 : t89 = t77*t13
388 2297800 : t93 = my_rho_1_3*my_rho
389 2297800 : t94 = 0.1e1_dp/t93
390 2297800 : t95 = t88*t94
391 : t98 = -0.240e3_dp*t48 - 0.216e3_dp*t5*d - 0.24e2_dp*t52&
392 : & *t5*t13*cf - 0.240e3_dp*t14*t48*cf - &
393 : 0.24e2_dp*t52*t2*t62 - 0.216e3_dp*t14*t5*cf &
394 : &*d + 0.10e2_dp/0.3e1_dp*t52*t70*t22 + 0.2e1_dp* &
395 : t14*t11*t22 + 0.10e2_dp/0.3e1_dp*t78*t80 + &
396 : 0.10e2_dp/0.3e1_dp*t14*t69*t22*d + 0.7e1_dp/ &
397 : 0.3e1_dp*t88*t89*t22 + 0.7e1_dp/0.3e1_dp*t95* &
398 2297800 : t80
399 2297800 : t99 = t98*t40
400 2297800 : t102 = 0.1e1_dp/t48
401 2297800 : t103 = a*t102
402 2297800 : t105 = 0.1e1_dp/t39/t38
403 2297800 : t106 = t37*t105
404 2297800 : t112 = my_rho_1_3*my_ndrho
405 : t123 = 0.6e1_dp*t14*t1*my_ndrho + 0.20e2_dp*t14*t112 &
406 : &*d + 0.14e2_dp*t14*t112*c + 0.14e2_dp*t14* &
407 2297800 : my_ndrho*c*d
408 2297800 : t124 = t123*t40
409 :
410 : e_rho(ii) = e_rho(ii) &
411 : - (0.5e1_dp/0.216e3_dp*t45*t41 - t4*t99/0.72e2_dp&
412 2297800 : & + t103*t106/0.108e3_dp)*sc
413 : e_ndrho(ii) = e_ndrho(ii) &
414 2297800 : + (t4*t124/0.72e2_dp)*sc
415 :
416 2297800 : IF (grad_deriv > 1) THEN
417 2297800 : t127 = 0.1e1_dp/t1/t6
418 2297800 : t128 = a*t127
419 2297800 : t133 = 0.1e1_dp/t7
420 2297800 : t134 = a*t133
421 2297800 : t161 = t3*t13
422 2297800 : t165 = 0.1e1_dp/t5
423 2297800 : t166 = t165*t13
424 2297800 : t173 = t52*t165
425 2297800 : t176 = t88*t102
426 2297800 : t184 = b*t87*c
427 2297800 : t185 = t102*t13
428 2297800 : t189 = t184*t44
429 : t192 = -0.560e3_dp*t93 - 0.432e3_dp*my_rho*d - 0.128e3_dp&
430 : & *t52*my_rho*t13*cf - 0.8e1_dp*t88*t1*t13* &
431 : cf - 0.560e3_dp*t14*t93*cf - 0.112e3_dp*t52*t1&
432 : & *t62 - 0.8e1_dp*t88*my_rho_1_3*t62 - 0.432e3_dp* &
433 : t14*my_rho*cf*d - 0.14e2_dp/0.9e1_dp*t52* &
434 : t161*t22 - 0.11e2_dp/0.9e1_dp*t88*t166*t22 - &
435 : 0.2e1_dp/0.3e1_dp*t14*t94*t22 - 0.20e2_dp/ &
436 : 0.9e1_dp*t173*t80 - 0.2e1_dp*t176*t80 - &
437 : 0.20e2_dp/0.9e1_dp*t14*t3*t22*d + 0.7e1_dp/ &
438 : 0.9e1_dp*t184*t185*t22 + 0.7e1_dp/0.9e1_dp* &
439 2297800 : t189*t80
440 2297800 : t193 = t192*t40
441 2297800 : t196 = t98*t105
442 2297800 : t199 = 0.1e1_dp/t6
443 2297800 : t200 = a*t199
444 2297800 : t201 = t39**2
445 2297800 : t202 = 0.1e1_dp/t201
446 2297800 : t203 = t37*t202
447 2297800 : t215 = t13*my_ndrho*d
448 : t227 = 0.20e2_dp/0.3e1_dp*t52*t70*my_ndrho + 0.4e1_dp* &
449 : t14*t11*my_ndrho + 0.20e2_dp/0.3e1_dp*t78*t215&
450 : & + 0.20e2_dp/0.3e1_dp*t14*t69*my_ndrho*d + &
451 : 0.14e2_dp/0.3e1_dp*t88*t89*my_ndrho + 0.14e2_dp &
452 2297800 : &/0.3e1_dp*t95*t215
453 2297800 : t228 = t227*t40
454 2297800 : t231 = t123*t105
455 : t245 = 0.6e1_dp*t14*t1 + 0.20e2_dp*t14*my_rho_1_3*d + &
456 : 0.14e2_dp*t14*my_rho_1_3*c + 0.14e2_dp*t14*c* &
457 2297800 : d
458 2297800 : t246 = t245*t40
459 :
460 : e_rho_rho(ii) = e_rho_rho(ii) &
461 : + (0.5e1_dp/0.81e2_dp*t128*t41 - 0.5e1_dp/0.108e3_dp&
462 : & *t45*t99 + t134*t106/0.27e2_dp + t4*t193/ &
463 : 0.72e2_dp - t103*t196/0.54e2_dp + t200*t203/ &
464 2297800 : 0.108e3_dp)*sc
465 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
466 : - (0.5e1_dp/0.216e3_dp*t45*t124 - t4*t228/ &
467 2297800 : 0.72e2_dp + t103*t231/0.108e3_dp)*sc
468 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
469 2297800 : + (t4*t246/0.72e2_dp)*sc
470 :
471 2297800 : IF (grad_deriv > 2) THEN
472 0 : t248 = t5**2
473 0 : t265 = 0.1e1_dp/t248
474 0 : t278 = t87**2
475 0 : t279 = b*t278
476 : t332 = -0.432e3_dp*d - 0.2240e4_dp/0.3e1_dp*my_rho_1_3 - &
477 : 0.74e2_dp/0.27e2_dp*t184*t127*t80 + 0.100e3_dp/ &
478 : 0.27e2_dp*t14*t44*t22*d + 0.7e1_dp/0.27e2_dp* &
479 : t279*t127*t13*t22 - 0.8e1_dp/0.3e1_dp*t184* &
480 : t70*cf - 0.2240e4_dp/0.3e1_dp*t14*my_rho_1_3* &
481 : cf - 0.48e2_dp*t88*t11*t13*cf - 0.944e3_dp/ &
482 : 0.3e1_dp*t52*t61 - 0.40e2_dp*t88*t69*t62 - &
483 : 0.656e3_dp/0.3e1_dp*t52*t11*t62 + 0.7e1_dp/ &
484 : 0.27e2_dp*t279*t265*t80 + 0.64e2_dp/0.27e2_dp* &
485 : t52*t44*t13*t22 - 0.8e1_dp/0.3e1_dp*t184*t77&
486 : & *t62 - 0.432e3_dp*t14*cf*d + 0.52e2_dp/ &
487 : 0.27e2_dp*t88*t199*t13*t22 - 0.20e2_dp/ &
488 : 0.9e1_dp*t184*t133*t13*t22 + 0.8e1_dp/0.9e1_dp&
489 : & *t14*t102*t22 + 0.100e3_dp/0.27e2_dp*t52*t199&
490 0 : & *t80 + 0.106e3_dp/0.27e2_dp*t88*t133*t80
491 :
492 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
493 : - (0.55e2_dp/0.243e3_dp*a/t1/t248*t41 - 0.5e1_dp &
494 : &/0.27e2_dp*t128*t99 + 0.40e2_dp/0.243e3_dp*a/ &
495 : my_rho_1_3/t248*t106 + 0.5e1_dp/0.72e2_dp*t45* &
496 : t193 - t134*t196/0.9e1_dp + 0.7e1_dp/0.108e3_dp* &
497 : a*t265*t203 - t4*t332*t40/0.72e2_dp + t103* &
498 : t192*t105/0.36e2_dp - t200*t98*t202/0.36e2_dp &
499 0 : &+ t128*t37/t201/t38/0.81e2_dp)*sc
500 : e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
501 : + (0.5e1_dp/0.81e2_dp*t128*t124 - 0.5e1_dp/ &
502 : 0.108e3_dp*t45*t228 + t134*t231/0.27e2_dp + t4*&
503 : & (-0.28e2_dp/0.9e1_dp*t52*t161*my_ndrho - &
504 : 0.22e2_dp/0.9e1_dp*t88*t166*my_ndrho - 0.4e1_dp &
505 : &/0.3e1_dp*t14*t94*my_ndrho - 0.40e2_dp/0.9e1_dp &
506 : &*t173*t215 - 0.4e1_dp*t176*t215 - 0.40e2_dp/ &
507 : 0.9e1_dp*t14*t3*my_ndrho*d + 0.14e2_dp/ &
508 : 0.9e1_dp*t184*t185*my_ndrho + 0.14e2_dp/0.9e1_dp&
509 : & *t189*t215)*t40/0.72e2_dp - t103*t227*t105/ &
510 0 : 0.54e2_dp + t200*t123*t202/0.108e3_dp)*sc
511 : e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
512 : - (0.5e1_dp/0.216e3_dp*t45*t246 - t4*(0.20e2_dp/ &
513 : 0.3e1_dp*t52*t70 + 0.4e1_dp*t14*t11 + 0.20e2_dp &
514 : &/0.3e1_dp*t52*t89*d + 0.20e2_dp/0.3e1_dp*t14* &
515 : t69*d + 0.14e2_dp/0.3e1_dp*t88*t89 + 0.14e2_dp/ &
516 : 0.3e1_dp*t88*t94*t13*d)*t40/0.72e2_dp + t103&
517 0 : & *t245*t105/0.108e3_dp)*sc
518 : END IF
519 : END IF
520 : END IF
521 : END IF
522 : END DO
523 :
524 : !$OMP END DO
525 :
526 : END SELECT
527 8277 : END SUBROUTINE lyp_lda_calc
528 :
529 : ! **************************************************************************************************
530 : !> \brief evaluates the becke 88 exchange functional for lsd
531 : !> \param rho_set ...
532 : !> \param deriv_set ...
533 : !> \param grad_deriv ...
534 : !> \param lyp_params input parameter (scaling)
535 : !> \par History
536 : !> 11.2003 created [fawzi]
537 : !> 01.2007 added scaling [Manuel Guidon]
538 : !> \author fawzi
539 : ! **************************************************************************************************
540 9120 : SUBROUTINE lyp_lsd_eval(rho_set, deriv_set, grad_deriv, lyp_params)
541 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
542 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
543 : INTEGER, INTENT(in) :: grad_deriv
544 : TYPE(section_vals_type), POINTER :: lyp_params
545 :
546 : CHARACTER(len=*), PARAMETER :: routineN = 'lyp_lsd_eval'
547 :
548 : INTEGER :: handle, npoints
549 : INTEGER, DIMENSION(2, 3) :: bo
550 : REAL(kind=dp) :: epsilon_drho, epsilon_rho, sc
551 3040 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
552 3040 : e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, &
553 3040 : e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, &
554 3040 : norm_drhob, rhoa, rhob
555 : TYPE(xc_derivative_type), POINTER :: deriv
556 :
557 3040 : CALL timeset(routineN, handle)
558 3040 : NULLIFY (deriv)
559 :
560 3040 : CALL section_vals_val_get(lyp_params, "scale_c", r_val=sc)
561 3040 : CALL cite_reference(Lee1988)
562 :
563 : CALL xc_rho_set_get(rho_set, &
564 : rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
565 : norm_drhob=norm_drhob, norm_drho=norm_drho, &
566 : rho_cutoff=epsilon_rho, &
567 3040 : drho_cutoff=epsilon_drho, local_bounds=bo)
568 3040 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
569 :
570 3040 : dummy => rhoa
571 3040 : e_0 => dummy
572 3040 : e_ra => dummy
573 3040 : e_rb => dummy
574 3040 : e_ndra_ra => dummy
575 3040 : e_ndra_rb => dummy
576 3040 : e_ndrb_ra => dummy
577 3040 : e_ndrb_rb => dummy
578 3040 : e_ndr_ndr => dummy
579 3040 : e_ndra_ndra => dummy
580 3040 : e_ndrb_ndrb => dummy
581 3040 : e_ndr => dummy
582 3040 : e_ndra => dummy
583 3040 : e_ndrb => dummy
584 3040 : e_ra_ra => dummy
585 3040 : e_ra_rb => dummy
586 3040 : e_rb_rb => dummy
587 3040 : e_ndr_ra => dummy
588 3040 : e_ndr_rb => dummy
589 :
590 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
591 3040 : allocate_deriv=.TRUE.)
592 3040 : CALL xc_derivative_get(deriv, deriv_data=e_0)
593 3040 : IF (grad_deriv >= 1) THEN
594 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
595 2956 : allocate_deriv=.TRUE.)
596 2956 : CALL xc_derivative_get(deriv, deriv_data=e_ra)
597 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
598 2956 : allocate_deriv=.TRUE.)
599 2956 : CALL xc_derivative_get(deriv, deriv_data=e_rb)
600 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
601 2956 : allocate_deriv=.TRUE.)
602 2956 : CALL xc_derivative_get(deriv, deriv_data=e_ndr)
603 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
604 2956 : allocate_deriv=.TRUE.)
605 2956 : CALL xc_derivative_get(deriv, deriv_data=e_ndra)
606 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
607 2956 : allocate_deriv=.TRUE.)
608 2956 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
609 : END IF
610 3040 : IF (grad_deriv >= 2) THEN
611 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
612 4 : allocate_deriv=.TRUE.)
613 4 : CALL xc_derivative_get(deriv, deriv_data=e_ra_ra)
614 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob], &
615 4 : allocate_deriv=.TRUE.)
616 4 : CALL xc_derivative_get(deriv, deriv_data=e_ra_rb)
617 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
618 4 : allocate_deriv=.TRUE.)
619 4 : CALL xc_derivative_get(deriv, deriv_data=e_rb_rb)
620 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhoa], &
621 4 : allocate_deriv=.TRUE.)
622 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndr_ra)
623 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rhob], &
624 4 : allocate_deriv=.TRUE.)
625 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndr_rb)
626 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
627 4 : allocate_deriv=.TRUE.)
628 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ra)
629 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhob], &
630 4 : allocate_deriv=.TRUE.)
631 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_rb)
632 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
633 4 : allocate_deriv=.TRUE.)
634 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb_rb)
635 : deriv => xc_dset_get_derivative(deriv_set, &
636 4 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
637 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndr_ndr)
638 : deriv => xc_dset_get_derivative(deriv_set, &
639 4 : [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
640 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndra_ndra)
641 : deriv => xc_dset_get_derivative(deriv_set, &
642 4 : [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
643 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb_ndrb)
644 : END IF
645 :
646 : !$OMP PARALLEL DEFAULT(NONE) &
647 : !$OMP SHARED(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob) &
648 : !$OMP SHARED(e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra) &
649 : !$OMP SHARED(e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb) &
650 : !$OMP SHARED(e_ndr, e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb) &
651 : !$OMP SHARED(e_ndr_ra, e_ndr_rb, grad_deriv) &
652 3040 : !$OMP SHARED(npoints, epsilon_rho, sc)
653 :
654 : CALL lyp_lsd_calc( &
655 : rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
656 : norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
657 : e_ndra_ra=e_ndra_ra, e_ndra_rb=e_ndra_rb, e_ndrb_ra&
658 : &=e_ndrb_ra, e_ndrb_rb=e_ndrb_rb, e_ndr_ndr=e_ndr_ndr, &
659 : e_ndra_ndra=e_ndra_ndra, e_ndrb_ndrb=e_ndrb_ndrb, e_ndr=e_ndr, &
660 : e_ndra=e_ndra, e_ndrb=e_ndrb, e_ra_ra=e_ra_ra, &
661 : e_ra_rb=e_ra_rb, e_rb_rb=e_rb_rb, e_ndr_ra=e_ndr_ra, &
662 : e_ndr_rb=e_ndr_rb, &
663 : grad_deriv=grad_deriv, npoints=npoints, &
664 : epsilon_rho=epsilon_rho, sc=sc)
665 :
666 : !$OMP END PARALLEL
667 :
668 3040 : CALL timestop(handle)
669 3040 : END SUBROUTINE lyp_lsd_eval
670 :
671 : ! **************************************************************************************************
672 : !> \brief low level calculation of the becke 88 exchange functional for lsd
673 : !> \param rhoa alpha spin density
674 : !> \param rhob beta spin density
675 : !> \param norm_drho || grad rhoa + grad rhob ||
676 : !> \param norm_drhoa || grad rhoa ||
677 : !> \param norm_drhob || grad rhob ||
678 : !> \param e_0 adds to it the local value of the functional
679 : !> \param e_ra e_*: derivative of the functional wrt. to the variables
680 : !> named where the * is.
681 : !> \param e_rb ...
682 : !> \param e_ndra_ra ...
683 : !> \param e_ndra_rb ...
684 : !> \param e_ndrb_ra ...
685 : !> \param e_ndrb_rb ...
686 : !> \param e_ndr_ndr ...
687 : !> \param e_ndra_ndra ...
688 : !> \param e_ndrb_ndrb ...
689 : !> \param e_ndr ...
690 : !> \param e_ndra ...
691 : !> \param e_ndrb ...
692 : !> \param e_ra_ra ...
693 : !> \param e_ra_rb ...
694 : !> \param e_rb_rb ...
695 : !> \param e_ndr_ra ...
696 : !> \param e_ndr_rb ...
697 : !> \param grad_deriv ...
698 : !> \param npoints ...
699 : !> \param epsilon_rho ...
700 : !> \param sc scaling parameter for correlation
701 : !> \par History
702 : !> 01.2004 created [fawzi]
703 : !> \author fawzi
704 : ! **************************************************************************************************
705 3040 : SUBROUTINE lyp_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
706 : e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, &
707 : e_ndra_ndra, e_ndrb_ndrb, e_ndr, &
708 : e_ndra, e_ndrb, e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb, &
709 : grad_deriv, npoints, epsilon_rho, sc)
710 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, &
711 : norm_drhob
712 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ndra_ra, e_ndra_rb, &
713 : e_ndrb_ra, e_ndrb_rb, e_ndr_ndr, e_ndra_ndra, e_ndrb_ndrb, e_ndr, e_ndra, e_ndrb, &
714 : e_ra_ra, e_ra_rb, e_rb_rb, e_ndr_ra, e_ndr_rb
715 : INTEGER, INTENT(in) :: grad_deriv, npoints
716 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sc
717 :
718 : INTEGER :: ii
719 : REAL(kind=dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rho, my_rhoa, my_rhob, t1, t100, &
720 : t101, t102, t103, t104, t107, t108, t109, t112, t115, t118, t12, t120, t124, t129, t13, &
721 : t130, t132, t135, t138, t14, t140, t142, t145, t15, t153, t155, t159, t16, t164, t168, &
722 : t17, t171, t176, t18, t181, t182, t185, t189, t194, t195, t198, t2, t20, t202, t205, &
723 : t206, t21, t210, t214, t215, t218, t22, t222, t228, t23, t232, t236, t238, t24, t243, &
724 : t249, t25, t252, t253, t255, t257, t26, t260, t265, t268, t27, t270, t29, t3, t30, t304, &
725 : t31, t310, t313, t316, t319, t322, t326, t329, t332, t334, t341, t35
726 : REAL(kind=dp) :: t354, t356, t360, t363, t37, t373, t376, t381, t39, t391, t4, t40, t408, &
727 : t41, t415, t419, t422, t430, t44, t445, t449, t45, t452, t46, t467, t47, t48, t483, t487, &
728 : t49, t490, t5, t50, t505, t515, t519, t52, t520, t54, t56, t57, t6, t61, t62, t64, t66, &
729 : t7, t72, t75, t78, t8, t82, t85, t86, t87, t88, t90, t92, t94, t95, t98, t99
730 :
731 3040 : cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
732 :
733 3040 : !$OMP DO
734 :
735 : DO ii = 1, npoints
736 54792107 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
737 54792107 : my_rhob = MAX(rhob(ii), 0.0_dp)
738 54792107 : my_rho = my_rhoa + my_rhob
739 54792107 : IF (my_rho > epsilon_rho) THEN
740 53280490 : my_ndrho = norm_drho(ii)
741 53280490 : my_ndrhoa = norm_drhoa(ii)
742 53280490 : my_ndrhob = norm_drhob(ii)
743 :
744 53280490 : t1 = my_rho**(0.1e1_dp/0.3e1_dp)
745 53280490 : t2 = 0.1e1_dp/t1
746 53280490 : t3 = d*t2
747 53280490 : t4 = 0.1e1_dp + t3
748 53280490 : t5 = 0.1e1_dp/t4
749 53280490 : t6 = a*t5
750 53280490 : t7 = my_rhoa*my_rhob
751 53280490 : t8 = 0.1e1_dp/my_rho
752 53280490 : t12 = a*b
753 53280490 : t13 = c*t2
754 53280490 : t14 = EXP(-t13)
755 53280490 : t15 = t12*t14
756 53280490 : t16 = my_rho**2
757 53280490 : t17 = t16*my_rho
758 53280490 : t18 = t1**2
759 53280490 : t20 = 0.1e1_dp/t18/t17
760 53280490 : t21 = t5*t20
761 53280490 : t22 = 2**(0.1e1_dp/0.3e1_dp)
762 53280490 : t23 = t22**2
763 53280490 : t24 = t23*cf
764 53280490 : t25 = my_rhoa**2
765 53280490 : t26 = my_rhoa**(0.1e1_dp/0.3e1_dp)
766 53280490 : t27 = t26**2
767 53280490 : t29 = my_rhob**2
768 53280490 : t30 = my_rhob**(0.1e1_dp/0.3e1_dp)
769 53280490 : t31 = t30**2
770 53280490 : t35 = 0.8e1_dp*t24*(t27*t25 + t31*t29)
771 53280490 : t37 = t3*t5
772 : t39 = 0.47e2_dp/0.18e2_dp - 0.7e1_dp/0.18e2_dp*t13 - &
773 53280490 : 0.7e1_dp/0.18e2_dp*t37
774 53280490 : t40 = my_ndrho**2
775 53280490 : t41 = t39*t40
776 53280490 : t44 = 0.5e1_dp/0.2e1_dp - t13/0.18e2_dp - t37/0.18e2_dp
777 53280490 : t45 = my_ndrhoa**2
778 53280490 : t46 = my_ndrhob**2
779 53280490 : t47 = t45 + t46
780 53280490 : t48 = t44*t47
781 53280490 : t49 = t13 + t37 - 0.11e2_dp
782 53280490 : t50 = my_rhoa*t8
783 53280490 : t52 = my_rhob*t8
784 53280490 : t54 = t50*t45 + t52*t46
785 53280490 : t56 = t49*t54/0.9e1_dp
786 53280490 : t57 = t35 + t41 - t48 - t56
787 53280490 : t61 = 0.2e1_dp/0.3e1_dp*t16
788 53280490 : t62 = t61 - t25
789 53280490 : t64 = t61 - t29
790 : t66 = t7*t57 - 0.2e1_dp/0.3e1_dp*t16*t40 + t62*t46 + &
791 53280490 : t64*t45
792 :
793 : IF (my_rho > epsilon_rho) THEN
794 : e_0(ii) = e_0(ii) &
795 53280490 : - (0.4e1_dp*t6*t7*t8 + t15*t21*t66)*sc
796 : END IF
797 : !--------
798 :
799 53280490 : IF (grad_deriv /= 0) THEN
800 47408114 : t72 = t27*my_rhoa
801 47408114 : t75 = t49*t8
802 47408114 : t78 = 0.64e2_dp/0.3e1_dp*t24*t72 - t75*t45/0.9e1_dp
803 47408114 : t82 = my_rhob*t57 + t7*t78 - 0.2e1_dp*my_rhoa*t46
804 47408114 : t85 = t4**2
805 47408114 : t86 = 0.1e1_dp/t85
806 47408114 : t87 = a*t86
807 47408114 : t88 = t87*my_rhoa
808 47408114 : t90 = 0.1e1_dp/t1/t16
809 47408114 : t92 = my_rhob*t90*d
810 47408114 : t94 = 0.4e1_dp/0.3e1_dp*t88*t92
811 47408114 : t95 = 0.1e1_dp/t16
812 47408114 : t98 = 0.4e1_dp*t6*t7*t95
813 47408114 : t99 = t12*c
814 47408114 : t100 = t16**2
815 47408114 : t101 = t100*my_rho
816 47408114 : t102 = 0.1e1_dp/t101
817 47408114 : t103 = t102*t14
818 47408114 : t104 = t5*t66
819 47408114 : t107 = t99*t103*t104/0.3e1_dp
820 47408114 : t108 = t86*t102
821 47408114 : t109 = t66*d
822 47408114 : t112 = t15*t108*t109/0.3e1_dp
823 47408114 : t115 = t5/t18/t100
824 47408114 : t118 = 0.11e2_dp/0.3e1_dp*t15*t115*t66
825 47408114 : t120 = 0.1e1_dp/t1/my_rho
826 47408114 : t124 = d**2
827 47408114 : t129 = c*t120 + d*t120*t5 - t124/t18/my_rho*t86
828 47408114 : t130 = 0.7e1_dp/0.54e2_dp*t129
829 47408114 : t132 = t129/0.54e2_dp
830 47408114 : t135 = -t129/0.3e1_dp
831 47408114 : t138 = my_rhoa*t95
832 47408114 : t140 = my_rhob*t95
833 47408114 : t142 = -t138*t45 - t140*t46
834 : t145 = t130*t40 - t132*t47 - t135*t54/0.9e1_dp - t49* &
835 47408114 : t142/0.9e1_dp
836 : t153 = t7*t145 - 0.4e1_dp/0.3e1_dp*my_rho*t40 + &
837 : 0.4e1_dp/0.3e1_dp*my_rho*t46 + 0.4e1_dp/0.3e1_dp&
838 47408114 : & *my_rho*t45
839 47408114 : t155 = t15*t21*t153
840 :
841 : e_ra(ii) = e_ra(ii) &
842 : - (0.4e1_dp*t6*t52 + t15*t21*t82 + t94 - t98 + &
843 47408114 : t107 + t112 - t118 + t155)*sc
844 :
845 47408114 : t159 = t31*my_rhob
846 47408114 : t164 = 0.64e2_dp/0.3e1_dp*t24*t159 - t75*t46/0.9e1_dp
847 47408114 : t168 = my_rhoa*t57 + t7*t164 - 0.2e1_dp*my_rhob*t45
848 :
849 : e_rb(ii) = e_rb(ii) &
850 : - (0.4e1_dp*t6*t50 + t15*t21*t168 + t94 - t98 + &
851 47408114 : t107 + t112 - t118 + t155)*sc
852 :
853 47408114 : t171 = t39*my_ndrho
854 47408114 : t176 = 0.2e1_dp*t7*t171 - 0.4e1_dp/0.3e1_dp*t16*my_ndrho
855 :
856 : e_ndr(ii) = e_ndr(ii) &
857 47408114 : - (t15*t21*t176)*sc
858 :
859 47408114 : t181 = t49*my_rhoa
860 47408114 : t182 = t8*my_ndrhoa
861 47408114 : t185 = -0.2e1_dp*t44*my_ndrhoa - 0.2e1_dp/0.9e1_dp*t181*t182
862 47408114 : t189 = t7*t185 + 0.2e1_dp*t64*my_ndrhoa
863 :
864 : e_ndra(ii) = e_ndra(ii) &
865 47408114 : - (t15*t21*t189)*sc
866 :
867 47408114 : t194 = t49*my_rhob
868 47408114 : t195 = t8*my_ndrhob
869 47408114 : t198 = -0.2e1_dp*t44*my_ndrhob - 0.2e1_dp/0.9e1_dp*t194*t195
870 47408114 : t202 = t7*t198 + 0.2e1_dp*t62*my_ndrhob
871 :
872 : e_ndrb(ii) = e_ndrb(ii) &
873 47408114 : - (t15*t21*t202)*sc
874 :
875 47408114 : IF (grad_deriv > 1) THEN
876 79625 : t205 = t100*t16
877 79625 : t206 = 0.1e1_dp/t205
878 79625 : t210 = 0.26e2_dp/0.9e1_dp*t99*t206*t14*t104
879 79625 : t214 = 0.2e1_dp/0.3e1_dp*t99*t103*t5*t153
880 79625 : t215 = c**2
881 79625 : t218 = 0.1e1_dp/t1/t205
882 79625 : t222 = t12*t215*t218*t14*t104/0.9e1_dp
883 79625 : t228 = 0.2e1_dp/0.9e1_dp*t12*c*t218*t14*t86*t109
884 79625 : t232 = 0.26e2_dp/0.9e1_dp*t15*t86*t206*t109
885 79625 : t236 = 0.2e1_dp/0.3e1_dp*t15*t108*t153*d
886 79625 : t238 = 0.1e1_dp/t85/t4
887 79625 : t243 = 0.2e1_dp/0.9e1_dp*t15*t238*t218*t66*t124
888 79625 : t249 = 0.154e3_dp/0.9e1_dp*t15*t5/t18/t101*t66
889 79625 : t252 = 0.22e2_dp/0.3e1_dp*t15*t115*t153
890 79625 : t253 = t6*t140
891 79625 : t255 = t87*t92
892 79625 : t257 = c*t90
893 79625 : t260 = d*t90*t5
894 79625 : t265 = t124/t18/t16*t86
895 79625 : t268 = 0.1e1_dp/t17
896 79625 : t270 = t124*d*t268*t238
897 : t304 = t15*t21*(t7*((-0.14e2_dp/0.81e2_dp*t257 - &
898 : 0.14e2_dp/0.81e2_dp*t260 + 0.7e1_dp/0.27e2_dp* &
899 : t265 - 0.7e1_dp/0.81e2_dp*t270)*t40 - (-0.2e1_dp/ &
900 : 0.81e2_dp*t257 - 0.2e1_dp/0.81e2_dp*t260 + t265/ &
901 : 0.27e2_dp - t270/0.81e2_dp)*t47 - (0.4e1_dp/ &
902 : 0.9e1_dp*t257 + 0.4e1_dp/0.9e1_dp*t260 - 0.2e1_dp &
903 : &/0.3e1_dp*t265 + 0.2e1_dp/0.9e1_dp*t270)*t54/ &
904 : 0.9e1_dp - 0.2e1_dp/0.9e1_dp*t135*t142 - t49*&
905 : & (0.2e1_dp*my_rhoa*t268*t45 + 0.2e1_dp*my_rhob* &
906 : t268*t46)/0.9e1_dp) - 0.4e1_dp/0.3e1_dp*t40 + &
907 79625 : 0.4e1_dp/0.3e1_dp*t46 + 0.4e1_dp/0.3e1_dp*t45)
908 79625 : t310 = 0.40e2_dp/0.9e1_dp*t88*my_rhob/t1/t17*d
909 79625 : t313 = my_rhob*t20
910 79625 : t316 = 0.8e1_dp/0.9e1_dp*a*t238*my_rhoa*t313*t124
911 79625 : t319 = 0.8e1_dp*t6*t7*t268
912 79625 : t322 = t99*t103*t5*t82
913 79625 : t326 = t15*t108*t82*d
914 79625 : t329 = t15*t115*t82
915 79625 : t332 = t135*t8
916 79625 : t334 = t49*t95
917 : t341 = t15*t21*(my_rhob*t145 + t7*(-t332*t45/ &
918 79625 : 0.9e1_dp + t334*t45/0.9e1_dp))
919 :
920 : e_ra_ra(ii) = e_ra_ra(ii) &
921 : + (t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 + &
922 : t252 + 0.8e1_dp*t253 - 0.8e1_dp/0.3e1_dp*t255 - &
923 : t304 + t310 - t316 - t319 - 0.2e1_dp/0.3e1_dp*t322 - &
924 : 0.2e1_dp/0.3e1_dp*t326 + 0.22e2_dp/0.3e1_dp*t329&
925 : & - 0.2e1_dp*t341 - t15*t21*(0.2e1_dp*my_rhob* &
926 : t78 + 0.320e3_dp/0.9e1_dp*t72*my_rhob*t24 - &
927 79625 : 0.2e1_dp*t46))*sc
928 :
929 79625 : t354 = t99*t103*t5*t168
930 79625 : t356 = t6*t138
931 79625 : t360 = t87*my_rhoa*t90*d
932 79625 : t363 = t15*t115*t168
933 : t373 = t15*t21*(my_rhoa*t145 + t7*(-t332*t46/ &
934 79625 : 0.9e1_dp + t334*t46/0.9e1_dp))
935 79625 : t376 = t15*t108*t168*d
936 :
937 : e_rb_rb(ii) = e_rb_rb(ii) &
938 : + (t210 - t214 - t222 - t228 + t232 - t236 - t243 - t249 + &
939 : t252 - t304 + t310 - t316 - t319 + 0.8e1_dp*t356 - &
940 : 0.8e1_dp/0.3e1_dp*t360 + 0.22e2_dp/0.3e1_dp* &
941 : t363 - 0.2e1_dp*t373 - 0.2e1_dp/0.3e1_dp*t354 - &
942 : 0.2e1_dp/0.3e1_dp*t376 - t15*t21*(0.2e1_dp* &
943 : my_rhoa*t164 + 0.320e3_dp/0.9e1_dp*my_rhoa* &
944 79625 : t159*t24 - 0.2e1_dp*t45))*sc
945 :
946 : t381 = -t354/0.3e1_dp + 0.4e1_dp*t356 - 0.4e1_dp/0.3e1_dp&
947 : & *t360 + 0.11e2_dp/0.3e1_dp*t363 - t373 - t341 - &
948 : t376/0.3e1_dp + t310 - 0.4e1_dp*t6*t8 + 0.4e1_dp* &
949 79625 : t253 + t210 - t214 - t222
950 : t391 = -t228 + t232 - t236 - t243 - t249 + t252 - 0.4e1_dp/ &
951 : 0.3e1_dp*t255 - t304 - t316 - t319 - t322/0.3e1_dp - &
952 : t326/0.3e1_dp + 0.11e2_dp/0.3e1_dp*t329 - t15* &
953 : t21*(t35 + t41 - t48 - t56 + my_rhob*t164 + my_rhoa &
954 79625 : &*t78)
955 :
956 : e_ra_rb(ii) = e_ra_rb(ii) &
957 79625 : + (t381 + t391)*sc
958 :
959 79625 : t408 = t12*t14*t5
960 79625 : t415 = t99*t103*t5*t176/0.3e1_dp
961 79625 : t419 = t15*t108*t176*d/0.3e1_dp
962 79625 : t422 = 0.11e2_dp/0.3e1_dp*t15*t115*t176
963 : t430 = t15*t21*(0.2e1_dp*t7*t130*my_ndrho - 0.8e1_dp &
964 79625 : &/0.3e1_dp*my_rho*my_ndrho)
965 :
966 : e_ndr_ra(ii) = e_ndr_ra(ii) &
967 : - (0.2e1_dp*t408*t313*t171 + t415 + t419 - t422 + &
968 79625 : t430)*sc
969 : e_ndr_rb(ii) = e_ndr_rb(ii) &
970 : - (0.2e1_dp*t408*t20*my_rhoa*t171 + t415 + t419 - &
971 79625 : t422 + t430)*sc
972 :
973 79625 : t445 = t99*t103*t5*t189/0.3e1_dp
974 79625 : t449 = t15*t108*t189*d/0.3e1_dp
975 79625 : t452 = 0.11e2_dp/0.3e1_dp*t15*t115*t189
976 : t467 = t15*t21*(t7*(-0.2e1_dp*t132*my_ndrhoa - &
977 : 0.2e1_dp/0.9e1_dp*t135*my_rhoa*t182 + 0.2e1_dp/ &
978 : 0.9e1_dp*t181*t95*my_ndrhoa) + 0.8e1_dp/0.3e1_dp&
979 79625 : & *my_rho*my_ndrhoa)
980 :
981 : e_ndra_ra(ii) = e_ndra_ra(ii) &
982 : - (t15*t21*(my_rhob*t185 - 0.2e1_dp/0.9e1_dp*t7&
983 79625 : & *t75*my_ndrhoa) + t445 + t449 - t452 + t467)*sc
984 : e_ndra_rb(ii) = e_ndra_rb(ii) &
985 : - (t15*t21*(my_rhoa*t185 - 0.4e1_dp*my_rhob* &
986 79625 : my_ndrhoa) + t445 + t449 - t452 + t467)*sc
987 :
988 79625 : t483 = t99*t103*t5*t202/0.3e1_dp
989 79625 : t487 = t15*t108*t202*d/0.3e1_dp
990 79625 : t490 = 0.11e2_dp/0.3e1_dp*t15*t115*t202
991 : t505 = t15*t21*(t7*(-0.2e1_dp*t132*my_ndrhob - &
992 : 0.2e1_dp/0.9e1_dp*t135*my_rhob*t195 + 0.2e1_dp/ &
993 : 0.9e1_dp*t194*t95*my_ndrhob) + 0.8e1_dp/0.3e1_dp&
994 79625 : & *my_rho*my_ndrhob)
995 :
996 : e_ndrb_ra(ii) = e_ndrb_ra(ii) &
997 : - (t15*t21*(my_rhob*t198 - 0.4e1_dp*my_rhoa* &
998 79625 : my_ndrhob) + t483 + t487 - t490 + t505)*sc
999 : e_ndrb_rb(ii) = e_ndrb_rb(ii) &
1000 : - (t15*t21*(my_rhoa*t198 - 0.2e1_dp/0.9e1_dp*t7&
1001 79625 : & *t75*my_ndrhob) + t483 + t487 - t490 + t505)*sc
1002 :
1003 79625 : t515 = 0.4e1_dp/0.3e1_dp*t16
1004 :
1005 : e_ndr_ndr(ii) = e_ndr_ndr(ii) &
1006 79625 : - (t15*t21*(0.2e1_dp*t7*t39 - t515))*sc
1007 :
1008 79625 : t519 = t13/0.9e1_dp
1009 79625 : t520 = t37/0.9e1_dp
1010 :
1011 : e_ndra_ndra(ii) = e_ndra_ndra(ii) &
1012 : - (t15*t21*(t7*(-0.5e1_dp + t519 + t520 - 0.2e1_dp &
1013 79625 : &/0.9e1_dp*t181*t8) + t515 - 0.2e1_dp*t29))*sc
1014 :
1015 : e_ndrb_ndrb(ii) = e_ndrb_ndrb(ii) &
1016 : - (t15*t21*(t7*(-0.5e1_dp + t519 + t520 - 0.2e1_dp &
1017 79625 : &/0.9e1_dp*t194*t8) + t515 - 0.2e1_dp*t25))*sc
1018 : END IF
1019 : END IF
1020 : END IF
1021 : END DO
1022 : !$OMP END DO
1023 :
1024 3040 : END SUBROUTINE lyp_lsd_calc
1025 :
1026 : END MODULE xc_lyp
|
