Skip to content

Commit b3ea54b

Browse files
simensrsimensr
authored
Merge pull request #87 from bast/doc
Refactor the Fortran example
2 parents a990c9c + 341ba22 commit b3ea54b

File tree

3 files changed

+338
-439
lines changed

3 files changed

+338
-439
lines changed

‎test/CMakeLists.txt‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -14,22 +14,22 @@ target_compile_options(testall
1414
add_test(testall ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/testall)
1515

1616
if(ENABLE_FC_SUPPORT)
17-
add_executable(kernel_example
18-
kernel_example.f90
17+
add_executable(fortran_example
18+
examples/fortran_example.f90
1919
${PROJECT_SOURCE_DIR}/api/xcfun.F90
2020
)
21-
target_link_libraries(kernel_example
21+
target_link_libraries(fortran_example
2222
XCFun
2323
)
24-
target_compile_options(kernel_example
24+
target_compile_options(fortran_example
2525
PRIVATE
2626
"$<$<CONFIG:Debug>:${XCFun_Fortran_FLAGS_DEBUG}>"
2727
"$<$<CONFIG:Release>:${XCFun_Fortran_FLAGS_RELEASE}>"
2828
"$<$<BOOL:${ENABLE_64BIT_INTEGERS}>:${XCFun_64BIT_INTEGERS_FLAGS}>"
2929
"$<$<BOOL:${ENABLE_CODE_COVERAGE}>:${XCFun_Fortran_FLAGS_COVERAGE}>"
3030
)
31-
set_target_properties(kernel_example PROPERTIES LINKER_LANGUAGE Fortran)
32-
add_test(kernel_example ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/kernel_example)
31+
set_target_properties(fortran_example PROPERTIES LINKER_LANGUAGE Fortran)
32+
add_test(fortran_example ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/fortran_example)
3333
endif()
3434

3535
if(ENABLE_PYTHON_INTERFACE)

‎test/examples/fortran_example.f90‎

Lines changed: 332 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,332 @@
1+
!
2+
! XCFun, an arbitrary order exchange-correlation library
3+
! Copyright (C) 2018 Ulf Ekström and contributors.
4+
!
5+
! This file is part of XCFun.
6+
!
7+
! This Source Code Form is subject to the terms of the Mozilla Public
8+
! License, v. 2.0. If a copy of the MPL was not distributed with this
9+
! file, You can obtain one at http://mozilla.org/MPL/2.0/.
10+
!
11+
! For information on the complete list of contributors to the
12+
! XCFun library, see: <https://xcfun.readthedocs.io/>
13+
program xc_example
14+
15+
! This example contains calls to f90 interface routines that are needed to "talk
16+
! to" the xcfun library and demonstrates how to use them.
17+
18+
! We will compute the kernel for an unpolarized system using total density and
19+
! the gradient components as the variables. These are linear in the density
20+
! matrix, which helps the code using the results from xcfun.
21+
22+
use xcfun, only: XC_CONTRACTED, &
23+
XC_N_NX_NY_NZ, &
24+
xcfun_splash, &
25+
xc_new_functional, &
26+
xc_set, &
27+
xc_eval_setup, &
28+
xc_eval, &
29+
xc_free_functional
30+
31+
implicit none
32+
33+
! we consider only one grid point
34+
integer, parameter :: num_grid_points = 1
35+
36+
! we will use XC_N_NX_NY_NZ
37+
! N: density
38+
! NX: x-gradient of the density
39+
! NY: y-gradient of the density
40+
! NZ: z-gradient of the density
41+
integer, parameter :: num_density_variables = 4
42+
43+
character(1000) :: text
44+
integer :: id, order, ierr, ipoint
45+
integer :: vector_length
46+
real(8) :: res
47+
48+
real(8), allocatable :: density(:, :, :)
49+
50+
51+
! print some info and copyright about the library
52+
! please always include this info in your code
53+
call xcfun_splash(text)
54+
print *, text(1:len_trim(text))
55+
56+
! create a new functional
57+
! we need this for interacting with the library
58+
id = xc_new_functional()
59+
60+
! in this example we use PBE
61+
print *, 'Setting up PBE'
62+
ierr = xc_set(id, 'pbe', 1.0d0)
63+
call assert(ierr == 0, "functional name not recognized")
64+
65+
66+
!-----------------------------------------------------------------------------
67+
! in the first example we compute the XC energy ("order 0 derivative")
68+
69+
order = 0
70+
! XC_CONTRACTED here has nothing to do with contracted basis sets
71+
! it means we will evaluated in the XC_CONTRACTED mode and
72+
! internally contract functional derivatives with the density taylor expansion
73+
! in other words: we will not have to explicitly assemble/contract partial
74+
! derivatives outside of XCFun
75+
ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
76+
call assert(ierr == 0, "xc_eval_setup failed")
77+
78+
vector_length = 2**order
79+
allocate(density(vector_length, num_density_variables, num_grid_points))
80+
density = 0.0d0
81+
do ipoint = 1, num_grid_points
82+
! we use fantasy values here
83+
density(1, 1, ipoint) = 1.0d0 ! n
84+
density(1, 2, ipoint) = 2.0d0 ! nabla_x n
85+
density(1, 3, ipoint) = 3.0d0 ! nabla_y n
86+
density(1, 4, ipoint) = 4.0d0 ! nabla_z n
87+
end do
88+
89+
res = derivative(id, &
90+
num_grid_points, &
91+
num_density_variables, &
92+
vector_length, &
93+
density)
94+
deallocate(density)
95+
print *, 'The XC energy density is', res
96+
97+
! compare with reference
98+
call assert((abs(-0.86494159400066051d0 - res) < 1.0d-6), &
99+
"derivatives do not match reference numbers")
100+
101+
102+
!-----------------------------------------------------------------------------
103+
! now we will compute the first derivatives ('potential')
104+
! and contract them with the first order densities
105+
106+
order = 1
107+
ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
108+
call assert(ierr == 0, "xc_eval_setup failed")
109+
110+
vector_length = 2**order
111+
allocate(density(vector_length, num_density_variables, num_grid_points))
112+
density = 0.0d0
113+
do ipoint = 1, num_grid_points
114+
! we use fantasy values here
115+
density(1, 1, ipoint) = 1.0d0 ! n zeroth order
116+
density(1, 2, ipoint) = 2.0d0 ! nabla_x n zeroth order
117+
density(1, 3, ipoint) = 3.0d0 ! nabla_y n zeroth order
118+
density(1, 4, ipoint) = 4.0d0 ! nabla_z n zeroth order
119+
density(2, 1, ipoint) = 5.0d0 ! n first order
120+
density(2, 2, ipoint) = 6.0d0 ! nabla_x n first order
121+
density(2, 3, ipoint) = 7.0d0 ! nabla_y n first order
122+
density(2, 4, ipoint) = 8.0d0 ! nabla_z n first order
123+
end do
124+
125+
res = derivative(id, &
126+
num_grid_points, &
127+
num_density_variables, &
128+
vector_length, &
129+
density)
130+
deallocate(density)
131+
132+
! compare with reference
133+
call assert((abs(-5.1509916226154067d0 - res) < 1.0d-6), &
134+
"derivatives do not match reference numbers")
135+
136+
137+
!-----------------------------------------------------------------------------
138+
! now we will compute a particular partial derivative
139+
! within order = 1
140+
! we do this with a trick: we set the perturbed density for
141+
! the density variable of interest to 1, and set other perturbed
142+
! densities to 0
143+
144+
allocate(density(vector_length, num_density_variables, num_grid_points))
145+
density = 0.0d0
146+
do ipoint = 1, num_grid_points
147+
! we use fantasy values here
148+
density(1, 1, ipoint) = 1.0d0 ! n zeroth order
149+
density(1, 2, ipoint) = 2.0d0 ! nabla_x n zeroth order
150+
density(1, 3, ipoint) = 3.0d0 ! nabla_y n zeroth order
151+
density(1, 4, ipoint) = 4.0d0 ! nabla_z n zeroth order
152+
density(2, 1, ipoint) = 0.0d0
153+
density(2, 2, ipoint) = 0.0d0
154+
density(2, 3, ipoint) = 1.0d0 ! we differentiate wrt this variable
155+
density(2, 4, ipoint) = 0.0d0
156+
end do
157+
158+
res = derivative(id, &
159+
num_grid_points, &
160+
num_density_variables, &
161+
vector_length, &
162+
density)
163+
deallocate(density)
164+
165+
! compare with reference
166+
call assert((abs(-1.3470456737102541d-2 - res) < 1.0d-6), &
167+
"derivatives do not match reference numbers")
168+
169+
170+
!-----------------------------------------------------------------------------
171+
! now we try 2nd order
172+
173+
order = 2
174+
ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
175+
call assert(ierr == 0, "xc_eval_setup failed")
176+
177+
vector_length = 2**order
178+
allocate(density(vector_length, num_density_variables, num_grid_points))
179+
density = 0.0d0
180+
do ipoint = 1, num_grid_points
181+
! we use fantasy values here
182+
density(1, 1, ipoint) = 1.0d0 ! zeroth order
183+
density(1, 2, ipoint) = 2.0d0 ! zeroth order
184+
density(1, 3, ipoint) = 3.0d0 ! zeroth order
185+
density(1, 4, ipoint) = 4.0d0 ! zeroth order
186+
density(2, 1, ipoint) = 5.0d0 ! first order
187+
density(2, 2, ipoint) = 6.0d0 ! first order
188+
density(2, 3, ipoint) = 7.0d0 ! first order
189+
density(2, 4, ipoint) = 8.0d0 ! first order
190+
density(3, 1, ipoint) = 5.0d0 ! first order
191+
density(3, 2, ipoint) = 6.0d0 ! first order
192+
density(3, 3, ipoint) = 7.0d0 ! first order
193+
density(3, 4, ipoint) = 8.0d0 ! first order
194+
density(4, 1, ipoint) = 0.0d0 ! second order
195+
density(4, 2, ipoint) = 0.0d0 ! second order
196+
density(4, 3, ipoint) = 0.0d0 ! second order
197+
density(4, 4, ipoint) = 0.0d0 ! second order
198+
end do
199+
200+
res = derivative(id, &
201+
num_grid_points, &
202+
num_density_variables, &
203+
vector_length, &
204+
density)
205+
deallocate(density)
206+
207+
! compare with reference
208+
call assert((abs(-9.4927931153398468d0 - res) < 1.0d-6), &
209+
"derivatives do not match reference numbers")
210+
211+
212+
!-----------------------------------------------------------------------------
213+
! now we try 3nd order, contracted with perturbed densities
214+
215+
order = 3
216+
ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order)
217+
call assert(ierr == 0, "xc_eval_setup failed")
218+
219+
vector_length = 2**order
220+
allocate(density(vector_length, num_density_variables, num_grid_points))
221+
density = 0.0d0
222+
do ipoint = 1, num_grid_points
223+
! we use fantasy values here
224+
density(1, 1, ipoint) = 1.0d0 ! zeroth order
225+
density(1, 2, ipoint) = 2.0d0 ! zeroth order
226+
density(1, 3, ipoint) = 3.0d0 ! zeroth order
227+
density(1, 4, ipoint) = 4.0d0 ! zeroth order
228+
density(2, 1, ipoint) = 5.0d0 ! first order (1)
229+
density(2, 2, ipoint) = 6.0d0 ! first order (1)
230+
density(2, 3, ipoint) = 7.0d0 ! first order (1)
231+
density(2, 4, ipoint) = 8.0d0 ! first order (1)
232+
density(3, 1, ipoint) = 9.0d0 ! first order (2)
233+
density(3, 2, ipoint) = 10.0d0 ! first order (2)
234+
density(3, 3, ipoint) = 11.0d0 ! first order (2)
235+
density(3, 4, ipoint) = 12.0d0 ! first order (2)
236+
density(4, 1, ipoint) = 5.0d0 ! second order (depending on (1) and (2))
237+
density(4, 2, ipoint) = 6.0d0 ! second order (depending on (1) and (2))
238+
density(4, 3, ipoint) = 7.0d0 ! second order (depending on (1) and (2))
239+
density(4, 4, ipoint) = 8.0d0 ! second order (depending on (1) and (2))
240+
density(5, 1, ipoint) = 9.0d0 ! first order (3)
241+
density(5, 2, ipoint) = 10.0d0 ! first order (3)
242+
density(5, 3, ipoint) = 11.0d0 ! first order (3)
243+
density(5, 4, ipoint) = 12.0d0 ! first order (3)
244+
density(6, 1, ipoint) = 5.0d0 ! second order (depending on (1) and (3))
245+
density(6, 2, ipoint) = 6.0d0 ! second order (depending on (1) and (3))
246+
density(6, 3, ipoint) = 7.0d0 ! second order (depending on (1) and (3))
247+
density(6, 4, ipoint) = 8.0d0 ! second order (depending on (1) and (3))
248+
density(7, 1, ipoint) = 9.0d0 ! second order (depending on (2) and (3))
249+
density(7, 2, ipoint) = 10.0d0 ! second order (depending on (2) and (3))
250+
density(7, 3, ipoint) = 11.0d0 ! second order (depending on (2) and (3))
251+
density(7, 4, ipoint) = 12.0d0 ! second order (depending on (2) and (3))
252+
density(8, 1, ipoint) = 0.0d0 ! third order (depending on (1), (2) and (3))
253+
density(8, 2, ipoint) = 0.0d0 ! third order (depending on (1), (2) and (3))
254+
density(8, 3, ipoint) = 0.0d0 ! third order (depending on (1), (2) and (3))
255+
density(8, 4, ipoint) = 0.0d0 ! third order (depending on (1), (2) and (3))
256+
end do
257+
258+
res = derivative(id, &
259+
num_grid_points, &
260+
num_density_variables, &
261+
vector_length, &
262+
density)
263+
deallocate(density)
264+
265+
! compare with reference
266+
call assert((abs(47.091223089835331d0 - res) < 1.0d-6), &
267+
"derivatives do not match reference numbers")
268+
269+
270+
!-----------------------------------------------------------------------------
271+
! we are done and can release the memory
272+
273+
call xc_free_functional(id)
274+
275+
print *, 'Kernel test has ended properly!'
276+
277+
contains
278+
279+
real(8) function derivative(id, &
280+
num_grid_points, &
281+
num_density_variables, &
282+
vector_length, &
283+
density)
284+
! computes the derivative and takes care of offsetting
285+
286+
integer, intent(in) :: id
287+
integer, intent(in) :: num_grid_points
288+
integer, intent(in) :: num_density_variables
289+
integer, intent(in) :: vector_length
290+
real(8), intent(in) :: density(vector_length, num_density_variables, num_grid_points)
291+
292+
real(8), allocatable :: input_array(:, :)
293+
real(8), allocatable :: output_array(:, :)
294+
295+
integer :: ipoint
296+
297+
allocate(input_array(num_density_variables*vector_length, num_grid_points))
298+
allocate(output_array(vector_length, num_grid_points))
299+
300+
! put the densities into the right places
301+
! along the input array
302+
do ipoint = 1, num_grid_points
303+
input_array(:, ipoint) = (/density(:, 1, ipoint), &
304+
density(:, 2, ipoint), &
305+
density(:, 3, ipoint), &
306+
density(:, 4, ipoint)/)
307+
end do
308+
309+
call xc_eval(id, num_grid_points, input_array, output_array)
310+
311+
! The output_array holds a Taylor series expansion
312+
! and we pick here one particular element out of this array.
313+
derivative = output_array(vector_length, 1)
314+
315+
deallocate(input_array)
316+
deallocate(output_array)
317+
end function
318+
319+
320+
subroutine assert(predicate, error_message)
321+
322+
logical, intent(in) :: predicate
323+
character(*), intent(in) :: error_message
324+
325+
if (.not. predicate) then
326+
print *, 'ERROR:', error_message
327+
stop 1
328+
endif
329+
330+
end subroutine
331+
332+
end program

0 commit comments

Comments
 (0)