Skip to content

Commit 3cc5e13

Browse files
robertodrrobertodr
authored
Merge pull request #111 from stigrj/cpp-example-2
Cpp example
2 parents 5d7b374 + 003d696 commit 3cc5e13

File tree

5 files changed

+367
-0
lines changed

5 files changed

+367
-0
lines changed

‎.ci/example.sh‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
#!/usr/bin/env bash
2+
3+
set -eu -o pipefail
4+
5+
example=$1
6+
lang="${example/_*/}"
7+
var=${lang}_COMPILER
8+
compiler=${var}
9+
echo "-- Running example $example"
10+
(
11+
cd examples/"$example"
12+
cmake -H. -Bbuild -DXCFun_DIR="$HOME"/Software/xcfun/share/cmake/XCFun -DCMAKE_"${lang}"_COMPILER="${!compiler}"
13+
cmake --build build
14+
cd build
15+
ctest
16+
)
17+
echo "-- Done with example $example"

‎.travis.yml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,9 @@ script:
126126
- ctest --output-on-failure --verbose --parallel 2
127127
# Check that Python bindings are usable from install prefix
128128
- env PYTHONPATH=$HOME/Software/xcfun/lib/python:$PYTHONPATH python -c 'import xcfun; print(xcfun.xcfun_splash())'
129+
- cd -
130+
# Build C++ example
131+
- bash .ci/example.sh CXX_host
129132

130133
after_success:
131134
- |

‎examples/CXX_host/CMakeLists.txt‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
cmake_minimum_required(VERSION 3.11)
2+
3+
project(CXX_host LANGUAGES CXX)
4+
5+
set(CMAKE_CXX_STANDARD 11)
6+
set(CMAKE_CXX_STANDARD_REQUIRED ON)
7+
set(CMAKE_CXX_EXTENSIONS OFF)
8+
9+
find_package(XCFun CONFIG QUIET)
10+
if(TARGET XCFun::xcfun)
11+
get_property(_loc TARGET XCFun::xcfun PROPERTY LOCATION)
12+
message(STATUS "Found XCFun: ${_loc} (found version ${XCFun_VERSION})")
13+
else()
14+
message(STATUS "Suitable XCFun could not be located. Fetching and building!")
15+
include(FetchContent)
16+
FetchContent_Declare(xcfun_sources
17+
QUIET
18+
URL
19+
https://github.com/dftlibs/xcfun/archive/v2.0.0a2.tar.gz
20+
)
21+
22+
FetchContent_GetProperties(xcfun_sources)
23+
24+
set(CMAKE_INSTALL_PREFIX ${CMAKE_BINARY_DIR}/external)
25+
set(ENABLE_TESTALL FALSE CACHE BOOL "")
26+
set(XCFUN_PYTHON_INTERFACE FALSE)
27+
28+
if(NOT xcfun_sources_POPULATED)
29+
FetchContent_Populate(xcfun_sources)
30+
31+
add_subdirectory(
32+
${xcfun_sources_SOURCE_DIR}
33+
${xcfun_sources_BINARY_DIR}
34+
)
35+
endif()
36+
endif()
37+
38+
enable_testing()
39+
40+
foreach(_src example)
41+
add_executable(${_src} ${_src}.cpp)
42+
target_link_libraries(${_src}
43+
PUBLIC
44+
XCFun::xcfun
45+
)
46+
add_test(
47+
NAME ${_src}
48+
COMMAND $<TARGET_FILE:${_src}>
49+
)
50+
endforeach()

‎examples/CXX_host/README.md‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
To fetch and build XCFun as a dependency
2+
----------------------------------------
3+
4+
$ mkdir build
5+
$ cd build
6+
$ cmake ..
7+
$ make
8+
$ make test
9+
10+
11+
To use an existing XCFun installation
12+
-------------------------------------
13+
14+
mkdir build
15+
cd build
16+
cmake .. -DXCFun_DIR=<path/to/xcfun>/share/cmake/XCFun
17+
make
18+
make test

‎examples/CXX_host/example.cpp‎

Lines changed: 279 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,279 @@
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+
*/
14+
15+
#include <cmath>
16+
#include <cstdlib>
17+
#include <iostream>
18+
19+
#include <XCFun/xcfun.h>
20+
21+
// This example contains calls to C++ interface routines that are needed to
22+
// "talk to" the xcfun library and demonstrates how to use them.
23+
24+
// We will compute the kernel for an unpolarized system using total density and
25+
// the gradient components as the variables. These are linear in the density
26+
// matrix, which helps the code using the results from xcfun.
27+
28+
29+
// we consider only one grid point
30+
const int num_grid_points = 1;
31+
32+
// we will use XC_N_NX_NY_NZ
33+
// N: density
34+
// NX: x-gradient of the density
35+
// NY: y-gradient of the density
36+
// NZ: z-gradient of the density
37+
const int num_density_variables = 4;
38+
39+
// forward declare function
40+
double derivative(xc_functional &id,
41+
int vector_length,
42+
double density[][num_density_variables][num_grid_points]);
43+
44+
int main(int, char **) {
45+
// print some info and copyright about the library
46+
// please always include this info in your code
47+
std::cout << xcfun_splash() << std::endl;
48+
49+
// create a new functional
50+
// we need this for interacting with the library
51+
auto id = xc_new_functional();
52+
53+
{
54+
// in this example we use PBE
55+
std::cout << "Setting up PBE" << std::endl;
56+
auto ierr = xc_set(id, "pbe", 1.0);
57+
if (ierr) std::cout << "functional name not recognized" << std::endl;
58+
}
59+
60+
//----------------------------------------------------------------------------
61+
// in the first example we compute the XC energy ("order 0 derivative")
62+
{
63+
const auto order = 0;
64+
// XC_CONTRACTED here has nothing to do with contracted basis sets
65+
// it means we will evaluated in the XC_CONTRACTED mode and internally
66+
// contract functional derivatives with the density taylor expansion
67+
// in other words: we will not have to explicitly assemble/contract partial
68+
// derivatives outside of XCFun
69+
auto ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order);
70+
if (ierr) std::cout << "xc_eval_setup failed" << std::endl;
71+
72+
auto vector_length = 1 << order; // bit shift to get power of two: 2**order
73+
double density[vector_length][num_density_variables][num_grid_points];
74+
75+
for (auto i = 0; i < num_grid_points; i++) {
76+
// we use fantasy values here
77+
density[0][0][i] = 1.0; // n
78+
density[0][1][i] = 2.0; // nabla_x n
79+
density[0][2][i] = 3.0; // nabla_y n
80+
density[0][3][i] = 4.0; // nabla_z n
81+
}
82+
83+
auto res = derivative(id, vector_length, density);
84+
std::cout << "The XC energy density is " << res << std::endl;
85+
86+
// compare with reference
87+
auto diff = std::abs(-0.86494159400066051 - res);
88+
if (diff > 1.0e-6)
89+
std::cout << "derivatives do not match reference numbers" << std::endl;
90+
}
91+
92+
//----------------------------------------------------------------------------
93+
// now we will compute the first derivatives ('potential')
94+
// and contract them with the first order densities
95+
{
96+
const auto order = 1;
97+
auto ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order);
98+
if (ierr) std::cout << "xc_eval_setup failed" << std::endl;
99+
100+
auto vector_length = 1 << order; // bit shift to get power of two: 2**order
101+
double density[vector_length][num_density_variables][num_grid_points];
102+
103+
for (auto i = 0; i < num_grid_points; i++) {
104+
// we use fantasy values here
105+
density[0][0][i] = 1.0; // n zeroth order
106+
density[0][1][i] = 2.0; // nabla_x n zeroth order
107+
density[0][2][i] = 3.0; // nabla_y n zeroth order
108+
density[0][3][i] = 4.0; // nabla_z n zeroth order
109+
density[1][0][i] = 5.0; // n first order
110+
density[1][1][i] = 6.0; // nabla_x n first order
111+
density[1][2][i] = 7.0; // nabla_y n first order
112+
density[1][3][i] = 8.0; // nabla_z n first order
113+
}
114+
115+
auto res = derivative(id, vector_length, density);
116+
117+
// compare with reference
118+
auto diff = std::abs(-5.1509916226154067 - res);
119+
if (diff > 1.0e-6)
120+
std::cout << "derivatives do not match reference numbers" << std::endl;
121+
}
122+
123+
//----------------------------------------------------------------------------
124+
// now we will compute a particular partial derivative
125+
// within order = 1
126+
// we do this with a trick: we set the perturbed density for
127+
// the density variable of interest to 1, and set other perturbed
128+
// densities to 0
129+
{
130+
const auto order = 1;
131+
auto ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order);
132+
if (ierr) std::cout << "xc_eval_setup failed" << std::endl;
133+
134+
const auto vector_length = 1 << order; // bit shift to get 2**order
135+
double density[vector_length][num_density_variables][num_grid_points];
136+
137+
for (auto i = 0; i < num_grid_points; i++) {
138+
// we use fantasy values here
139+
density[0][0][i] = 1.0; // n zeroth order
140+
density[0][1][i] = 2.0; // nabla_x n zeroth order
141+
density[0][2][i] = 3.0; // nabla_y n zeroth order
142+
density[0][3][i] = 4.0; // nabla_z n zeroth order
143+
density[1][0][i] = 0.0;
144+
density[1][1][i] = 0.0;
145+
density[1][2][i] = 1.0; // we differentiate wrt this variable
146+
density[1][3][i] = 0.0;
147+
}
148+
149+
auto res = derivative(id, vector_length, density);
150+
151+
// compare with reference
152+
auto diff = std::abs(-0.013470456737102541 - res);
153+
if (diff > 1.0e-6)
154+
std::cout << "derivatives do not match reference numbers" << std::endl;
155+
}
156+
157+
//----------------------------------------------------------------------------
158+
// now we try 2nd order
159+
{
160+
const auto order = 2;
161+
auto ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order);
162+
if (ierr) std::cout << "xc_eval_setup failed" << std::endl;
163+
164+
const auto vector_length = 1 << order; // bit shift to get 2**order
165+
double density[vector_length][num_density_variables][num_grid_points];
166+
167+
for (auto i = 0; i < num_grid_points; i++) {
168+
// we use fantasy values here
169+
density[0][0][i] = 1.0; // zeroth order
170+
density[0][1][i] = 2.0; // zeroth order
171+
density[0][2][i] = 3.0; // zeroth order
172+
density[0][3][i] = 4.0; // zeroth order
173+
density[1][0][i] = 5.0; // first order
174+
density[1][1][i] = 6.0; // first order
175+
density[1][2][i] = 7.0; // first order
176+
density[1][3][i] = 8.0; // first order
177+
density[2][0][i] = 5.0; // first order
178+
density[2][1][i] = 6.0; // first order
179+
density[2][2][i] = 7.0; // first order
180+
density[2][3][i] = 8.0; // first order
181+
density[3][0][i] = 0.0; // second order
182+
density[3][1][i] = 0.0; // second order
183+
density[3][2][i] = 0.0; // second order
184+
density[3][3][i] = 0.0; // second order
185+
}
186+
187+
auto res = derivative(id, vector_length, density);
188+
189+
// compare with reference
190+
auto diff = std::abs(-9.4927931153398468 - res);
191+
if (diff > 1.0e-6)
192+
std::cout << "derivatives do not match reference numbers" << std::endl;
193+
}
194+
195+
//----------------------------------------------------------------------------
196+
// now we try 3nd order, contracted with perturbed densities
197+
{
198+
const auto order = 3;
199+
auto ierr = xc_eval_setup(id, XC_N_NX_NY_NZ, XC_CONTRACTED, order);
200+
if (ierr) std::cout << "xc_eval_setup failed" << std::endl;
201+
202+
auto vector_length = 1 << order; // bit shift to get power of two: 2**order
203+
double density[vector_length][num_density_variables][num_grid_points];
204+
205+
for (auto i = 0; i < num_grid_points; i++) {
206+
// we use fantasy values here
207+
density[0][0][i] = 1.0; // zeroth order
208+
density[0][1][i] = 2.0; // zeroth order
209+
density[0][2][i] = 3.0; // zeroth order
210+
density[0][3][i] = 4.0; // zeroth order
211+
density[1][0][i] = 5.0; // first order (1)
212+
density[1][1][i] = 6.0; // first order (1)
213+
density[1][2][i] = 7.0; // first order (1)
214+
density[1][3][i] = 8.0; // first order (1)
215+
density[2][0][i] = 9.0; // first order (2)
216+
density[2][1][i] = 10.0; // first order (2)
217+
density[2][2][i] = 11.0; // first order (2)
218+
density[2][3][i] = 12.0; // first order (2)
219+
density[3][0][i] = 5.0; // second order (depending on (1) and (2))
220+
density[3][1][i] = 6.0; // second order (depending on (1) and (2))
221+
density[3][2][i] = 7.0; // second order (depending on (1) and (2))
222+
density[3][3][i] = 8.0; // second order (depending on (1) and (2))
223+
density[4][0][i] = 9.0; // first order (3)
224+
density[4][1][i] = 10.0; // first order (3)
225+
density[4][2][i] = 11.0; // first order (3)
226+
density[4][3][i] = 12.0; // first order (3)
227+
density[5][0][i] = 5.0; // second order (depending on (1) and (3))
228+
density[5][1][i] = 6.0; // second order (depending on (1) and (3))
229+
density[5][2][i] = 7.0; // second order (depending on (1) and (3))
230+
density[5][3][i] = 8.0; // second order (depending on (1) and (3))
231+
density[6][0][i] = 9.0; // second order (depending on (2) and (3))
232+
density[6][1][i] = 10.0; // second order (depending on (2) and (3))
233+
density[6][2][i] = 11.0; // second order (depending on (2) and (3))
234+
density[6][3][i] = 12.0; // second order (depending on (2) and (3))
235+
density[7][0][i] = 0.0; // third order (depending on (1), (2) and (3))
236+
density[7][1][i] = 0.0; // third order (depending on (1), (2) and (3))
237+
density[7][2][i] = 0.0; // third order (depending on (1), (2) and (3))
238+
density[7][3][i] = 0.0; // third order (depending on (1), (2) and (3))
239+
}
240+
241+
auto res = derivative(id, vector_length, density);
242+
243+
// compare with reference
244+
auto diff = std::abs(47.091223089835331 - res);
245+
if (diff > 1.0e-6)
246+
std::cout << "derivatives do not match reference numbers" << std::endl;
247+
}
248+
249+
//----------------------------------------------------------------------------
250+
// we are done and can release the memory
251+
xc_free_functional(id);
252+
std::cout << "Kernel test has ended properly!" << std::endl;
253+
return EXIT_SUCCESS;
254+
}
255+
256+
// computes the derivative and takes care of offsetting
257+
double derivative(xc_functional &id,
258+
int vector_length,
259+
double density[][num_density_variables][num_grid_points]) {
260+
double inp_array[vector_length * num_density_variables];
261+
double out_array[num_grid_points][vector_length];
262+
263+
// put the densities into the right places
264+
// along the input array
265+
for (auto i = 0; i < num_grid_points; i++) {
266+
auto n = 0;
267+
for (auto j = 0; j < num_density_variables; j++) {
268+
for (auto k = 0; k < vector_length; k++) {
269+
inp_array[n++] = density[k][j][i];
270+
}
271+
}
272+
xc_eval(id, inp_array, out_array[i]);
273+
}
274+
275+
// The output_array holds a Taylor series expansion
276+
// and we pick here one particular element out of this array.
277+
return out_array[0][vector_length - 1];
278+
}
279+

0 commit comments

Comments
 (0)