From 911f12778949e0f39394e9f52b53f65eb8ac4df1 Mon Sep 17 00:00:00 2001 From: "Oleg E. Peil" Date: Wed, 9 Mar 2016 18:47:30 +0100 Subject: [PATCH] Ported C-style 'dos_tetra' to C++ In order to wrap the ATM routines by Python using TRIQS wrapping tools it is necessary to modify the interface to 'dos_tetra3d'. The major changes involved replacing direct NumPy arrays with TRIQS arrays which can be converted to Python arrays using library tools. Also, some small changes were necessary to port the functions from C99 complex numbers to C++ style. CMakeList is added to automatize building of the ATM library. --- c++/plovasp/atm/CMakeLists.txt | 9 + c++/plovasp/atm/{argsort.c => argsort.cpp} | 2 +- c++/plovasp/atm/{argsort.h => argsort.hpp} | 0 .../atm/{dos_tetra3d.c => dos_tetra3d.cpp} | 317 ++++++++---------- c++/plovasp/atm/dos_tetra3d.h | 42 --- c++/plovasp/atm/dos_tetra3d.hpp | 46 +++ 6 files changed, 200 insertions(+), 216 deletions(-) create mode 100644 c++/plovasp/atm/CMakeLists.txt rename c++/plovasp/atm/{argsort.c => argsort.cpp} (98%) rename c++/plovasp/atm/{argsort.h => argsort.hpp} (100%) rename c++/plovasp/atm/{dos_tetra3d.c => dos_tetra3d.cpp} (59%) delete mode 100644 c++/plovasp/atm/dos_tetra3d.h create mode 100644 c++/plovasp/atm/dos_tetra3d.hpp diff --git a/c++/plovasp/atm/CMakeLists.txt b/c++/plovasp/atm/CMakeLists.txt new file mode 100644 index 00000000..2f222969 --- /dev/null +++ b/c++/plovasp/atm/CMakeLists.txt @@ -0,0 +1,9 @@ +# Linking and include info +add_library(atm_c dos_tetra3d.hpp dos_tetra3d.cpp argsort.hpp argsort.cpp) +set_target_properties(atm_c PROPERTIES LINKER_LANGUAGE CXX) +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/c++/plovasp/atm ${TRIQS_INCLUDE_ALL}) + +#add_executable(test_atm test2py.cpp) +#target_link_libraries(test_atm atm_c) + +add_subdirectory(test) diff --git a/c++/plovasp/atm/argsort.c b/c++/plovasp/atm/argsort.cpp similarity index 98% rename from c++/plovasp/atm/argsort.c rename to c++/plovasp/atm/argsort.cpp index fad642d9..92fb86c2 100644 --- a/c++/plovasp/atm/argsort.c +++ b/c++/plovasp/atm/argsort.cpp @@ -18,7 +18,7 @@ * TRIQS. If not, see . * *******************************************************************************/ -#include +#include int cmp(const void *a, const void *b) { diff --git a/c++/plovasp/atm/argsort.h b/c++/plovasp/atm/argsort.hpp similarity index 100% rename from c++/plovasp/atm/argsort.h rename to c++/plovasp/atm/argsort.hpp diff --git a/c++/plovasp/atm/dos_tetra3d.c b/c++/plovasp/atm/dos_tetra3d.cpp similarity index 59% rename from c++/plovasp/atm/dos_tetra3d.c rename to c++/plovasp/atm/dos_tetra3d.cpp index 19893296..f688484b 100644 --- a/c++/plovasp/atm/dos_tetra3d.c +++ b/c++/plovasp/atm/dos_tetra3d.cpp @@ -18,15 +18,18 @@ * TRIQS. If not, see . * *******************************************************************************/ -#include +#include +#include +#include -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include "argsort.hpp" +#include "dos_tetra3d.hpp" -#include -#include -//#include "tetra.h" -#include "argsort.h" -#include "dos_tetra3d.h" +//#define __TETRA_DEBUG +#define __TETRA_ARRAY_VIEW + +using triqs::arrays::array; +using triqs::arrays::array_view; /*************************************************** @@ -35,6 +38,23 @@ Lambin et al., PRB 29, 6, 3430 (1984). ***************************************************/ +/// Main function +//#ifdef __TETRA_ARRAY_VIEW +//void tet_dos3d(double en, array_view& eigk, +// array_view& itt, int ntet, +// array& cti); +//#else +//void tet_dos3d(double en, array& eigk, +// array& itt, int ntet, +// array& cti); +//#endif + + +/// Internal functions +int dos_corner_weights(double en, double *eigs, int *inds, double *ci); +int dos_tet_weights(double en, double *eigs, int *inds, double *ct); +int dos_reorder(double en, double *e, int *inds); + static double F(double en, double e1, double e2, double e3, double e4); static double K2(double en, double e1, double e2, double e3); static double K1(double en, double e1, double e2); @@ -43,180 +63,141 @@ static void fun_dos_case1(double en, double *eigs, double *ci); static void fun_dos_case2(double en, double *eigs, double *ci); static void fun_dos_case3(double en, double *eigs, double *ci); -int dos_corner_weights(double en, double *eigs, int *inds, double *ci); -int dos_tet_weights(double en, double *eigs, int *inds, double *ct); +static const int NUM_TET_CORNERS = 4; +static const std::complex I(0.0, 1.0); +static const double small = 2.5e-2, tol = 1e-8; -const static int NUM_TET_CORNERS = 4; - -static PyMethodDef c_tetraMethods[] = { - {"dos_weights_3d", tetra_DOS3D, METH_VARARGS, - "C-implementation of the tetrahedron method for calculating DOS"}, - {NULL, NULL, 0, NULL} -}; - -PyMODINIT_FUNC -initc_atm_dos(void) +/* + Returns corner contributions to the DOS of a band +*/ +#ifdef __TETRA_ARRAY_VIEW +array_view dos_tetra_weights_3d(array_view eigk, double en, array_view itt) +#else +array dos_tetra_weights_3d(array eigk, double en, array itt) +#endif { - (void) Py_InitModule("c_atm_dos", c_tetraMethods); - import_array(); -} - -// -// Integration_Weights -// -static PyObject * -tetra_DOS3D(PyObject *self, PyObject *args) -{ - - PyArrayObject *py_eigk, *py_itt; // Input Numpy arrays - PyArrayObject *py_cti; // Output Numpy array -// npy_int64 *itt; // C-pointer to the 'itt' array - npy_intp *dims; // Dimensions of the output array 'rti' - double *cti; // C-pointer to 'eigk' and 'rti' arrays - double e; -// double et[4]; // Real part of eigenvalues -// int strd_itt[2], strd_eigk, strd_cti[2], it_cti; // Strides - int ntet; // Indices of sorted 'eigs'; number of tetrahedra + int ntet; /// Number of tetrahedra // Auxiliary variables and loop indices - int nd; -// Data-preparation part - if (!PyArg_ParseTuple(args, "O!dO!", - &PyArray_Type, &py_eigk, &e, &PyArray_Type, &py_itt)) - return NULL; - -// Sanity tests (the types are assumed to be checked in the python wrapper) -// nd = py_eigk->nd; - nd = PyArray_NDIM(py_eigk); - if (nd != 1) + if (first_dim(itt) != NUM_TET_CORNERS + 1) { - PyErr_SetString(PyExc_ValueError, " Array 'eigk' must be 1D"); - return NULL; + TRIQS_RUNTIME_ERROR << " The first dimension of 'itt' must be equal to 5"; } -// nd = py_itt->nd; - nd = PyArray_NDIM(py_itt); - if (nd != 2) - { - PyErr_SetString(PyExc_ValueError, " Array 'itt' must be 2D"); - return NULL; - } + ntet = second_dim(itt); -// nd = py_itt->dimensions[0]; - dims = PyArray_DIMS(py_itt); - if (dims[0] != 5) - { - PyErr_SetString(PyExc_ValueError, - " The first dimension of 'itt' must be equal to 5"); - return NULL; - } + array cti(NUM_TET_CORNERS, ntet); // Corner weights to be returned -// ntet = (int) py_itt->dimensions[1]; - ntet = dims[1]; - -// eigk = (double *)py_eigk->data; -// strd_eigk = py_eigk->strides[0] / sizeof(npy_float64); - -// itt = (npy_int64 *)py_itt->data; -// strd_itt[0] = py_itt->strides[0] / sizeof(npy_int64); -// strd_itt[1] = py_itt->strides[1] / sizeof(npy_int64); - -// Resulting array (the question is whether dims is copied or not?) - cti = (double *)malloc(NUM_TET_CORNERS * ntet * sizeof(double)); - dims = (npy_intp *)malloc(2 * sizeof(npy_intp)); - - dims[0] = NUM_TET_CORNERS; - dims[1] = ntet; - - py_cti = (PyArrayObject *)PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, cti); - -// strd_cti[0] = py_cti->strides[0] / sizeof(double); -// strd_cti[1] = py_cti->strides[1] / sizeof(double); - -// Main part: now we can fill the 'cti' array and it will be returned -// by 'py_cti' as a numpy array +// tet_dos3d(e, eigk, itt, ntet, cti); // -// Main function +// Main algorithm (transferred from 'tet_dos3d()') // -// tet_dos3d(e, eigk, strd_eigk, itt, ntet, strd_itt, cti, strd_cti); - tet_dos3d(e, py_eigk, py_itt, ntet, py_cti); - - return PyArray_Return(py_cti); -} - -//void tet_dos3d(double en, double *eigk, int strd_eigk, -// npy_int64 *itt, int ntet, int *strd_itt, -// double *cti, int *strd_cti) -void tet_dos3d(double en, PyArrayObject *py_eigk, - PyArrayObject *py_itt, int ntet, - PyArrayObject *py_cti) -{ double eigs[4], ci[4]; + int i, it, ik, inds[4], flag; -// **** DEBUG +#ifdef __TETRA_DEBUG double ct, ci_sum; +#endif // Loop over tetrahedra (triangles) for (it = 0; it < ntet; it++) { -// it_cti = it * strd_cti[1]; -// Sort eigenvalues and obtain indices of the sorted array -// eigs: sorted eigenvalues -// inds: index map for (i = 1; i < 5; i++) { -// ik = itt[i * strd_itt[0] + it * strd_itt[1]]; - ik = ((int *)PyArray_GETPTR2(py_itt, i, it))[0]; -// eigs[i - 1] = eigk[ik * strd_eigk]; - eigs[i - 1] = ((double *)PyArray_GETPTR1(py_eigk, ik))[0]; + ik = itt(i, it); + eigs[i - 1] = eigk(ik); } -// corner weights for a single tetrahedron +// Corner weights for a single tetrahedron dos_corner_weights(en, eigs, inds, ci); + +#ifdef __TETRA_DEBUG for(i = 0, ci_sum = 0.0; i < 4; i++) ci_sum += ci[i]; flag = dos_tet_weights(en, eigs, inds, &ct); if(fabs(ct - ci_sum) > tol) { - printf(" *** Error in weights: it = %d, flag = %d, en = %lf", it, flag, en); + std::cout << " *** Error in weights: it = " << it <<" flag = " << flag << ", en = " << en; for(i = 0; i < 4; i++) - printf(", e[%d] = %lf", i, eigs[i]); - printf(", c_diff = %le\n", fabs(ct - ci_sum)); - return; + std::cout << ", e[" << i << "] = " << eigs[i]; + std::cout << ", c_diff = " << fabs(ct - ci_sum) << std::endl; + + TRIQS_RUNTIME_ERROR << " Failed consistency check"; } -// printf(" it = %d, flag = %d", it, j); -// for(i = 0; i < 4; i++) -// printf(", e[%d] = %lf", i, eigs[i]); -// printf(", c_diff = %le\n", fabs(ct - ci_sum)); +#endif -// if(j < 4) -// { -// printf(" flag = %d, e = %lf", j, en); -// for(i = 0; i < 4; i++) -// printf(", eigs[%d] = %lf", i, eigs[i]); -// printf("\n"); -// printf(" ci = "); -// for(i = 0; i < 4; i++) -// printf(", %lf", ci[i]); -// printf("\n"); -// } - - for(i = 0; i < 4; i++) + for(i = 0; i < 4; i++) { -// j = inds[i] * strd_cti[0] + it_cti; -// cti[j] = ci[i]; - ((double *)PyArray_GETPTR2(py_cti, inds[i], it))[0] = ci[i]; + cti(inds[i], it) = ci[i]; } } // it = 1, ntet + + return array_view(cti); } + +//#ifdef __TETRA_ARRAY_VIEW +//void tet_dos3d(double en, array_view& eigk, +// array_view& itt, int ntet, +// array& cti) +//#else +//void tet_dos3d(double en, array& eigk, +// array& itt, int ntet, +// array& cti) +//#endif +//{ +// double eigs[4], ci[4]; +// +// int i, it, ik, inds[4], flag; +//#ifdef __TETRA_DEBUG +// double ct, ci_sum; +//#endif +// +//// Loop over tetrahedra (triangles) +// for (it = 0; it < ntet; it++) +// { +// for (i = 1; i < 5; i++) +// { +// ik = itt(i, it); +// eigs[i - 1] = eigk(ik); +// } +// +//// Corner weights for a single tetrahedron +// dos_corner_weights(en, eigs, inds, ci); +// +//#ifdef __TETRA_DEBUG +// for(i = 0, ci_sum = 0.0; i < 4; i++) +// ci_sum += ci[i]; +// +// flag = dos_tet_weights(en, eigs, inds, &ct); +// if(fabs(ct - ci_sum) > tol) +// { +// std::cout << " *** Error in weights: it = " << it <<" flag = " << flag << ", en = " << en; +// for(i = 0; i < 4; i++) +// std::cout << ", e[" << i << "] = " << eigs[i]; +// std::cout << ", c_diff = " << fabs(ct - ci_sum) << std::endl; +// return; +// } +//#endif +// +// for(i = 0; i < 4; i++) +// { +// cti(inds[i], it) = ci[i]; +// } +// +// } // it = 1, ntet +//} +/// Corner contributions to DOS int dos_corner_weights(double en, double *eigs, int *inds, double *ci) { int flag, i; +// Sort eigenvalues and obtain indices of the sorted array +// eigs: sorted eigenvalues +// inds: index map flag = dos_reorder(en, eigs, inds); switch(flag) @@ -251,11 +232,14 @@ int dos_corner_weights(double en, double *eigs, int *inds, return flag; } +/// Total (tetrahedron) contribution to DOS. +/// Here, it is calculated directly using an analytical formula. +/// This is mainly needed for debugging. int dos_tet_weights(double en, double *eigs, int *inds, double *ct) { double e1, e2, e3, e4; - double complex s; + std::complex s; int flag; flag = dos_reorder(en, eigs, inds); @@ -275,9 +259,9 @@ int dos_tet_weights(double en, double *eigs, int *inds, s = fmin(fabs(e1 - e2), fabs(e3 - e1)); s = fmin(fabs(s), fabs(e4 - e1)); s /= 100.0; - s = fmax(s, 1.0e-20) * I; + s = fmax(std::abs(s), 1.0e-20) * I; - *ct = 3.0 * creal((en - e1 + s) * (en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s) * (e4 - e1 + s))); + *ct = 3.0 * std::real((en - e1 + s) * (en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s) * (e4 - e1 + s))); } break; @@ -294,9 +278,9 @@ int dos_tet_weights(double en, double *eigs, int *inds, s = fmin(fabs(s), fabs(e4 - e1)); s = fmin(fabs(s), fabs(e4 - e2)); s /= 100.0; - s = fmax(s, 1.0e-20) * I; + s = fmax(std::abs(s), 1.0e-20) * I; - *ct = 3.0 * creal(( + *ct = 3.0 * std::real(( (e3 - en + s) * (en - e2 + s) / ((e4 - e2 + s) * (e3 - e2 + s) * (e3 - e1 + s)) + (e4 - en + s) * (en - e1 + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e3 - e1 + s)))); } @@ -311,9 +295,9 @@ int dos_tet_weights(double en, double *eigs, int *inds, s = fmin(fabs(e4 - e2), fabs(e4 - e1)); s = fmin(fabs(s), fabs(e4 - e3)); s /= 100.0; - s = fmax(s, 1.0e-20) * I; + s = fmax(std::abs(s), 1.0e-20) * I; - *ct = 3.0 * creal((e4 - en + s) * (e4 - en + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e4 - e3 + s))); + *ct = 3.0 * std::real((e4 - en + s) * (e4 - en + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e4 - e3 + s))); } break; @@ -333,6 +317,8 @@ int dos_tet_weights(double en, double *eigs, int *inds, return flag; } +/// Sorts eigenvalues and also determines eigenvalue degeneracies. +/// Returns a case number corresponding to a combination of degeneracies. int dos_reorder(double en, double *e, int *inds) { double *ptrs[4], e_tmp[4]; @@ -358,13 +344,6 @@ int dos_reorder(double en, double *e, int *inds) static void fun_dos_case1(double en, double *eigs, double *ci) { double e1, e2, e3, e4; -// int i; - -// if(fabs(eigs[1] - eigs[0]) < tol) -// { -// for(i = 0; i < 4; i++) ci[i] = 0.0; -// return; -// } e1 = eigs[0]; e2 = eigs[1]; @@ -386,12 +365,6 @@ static void fun_dos_case2(double en, double *eigs, double *ci) { double e1, e2, e3, e4; -// if(fabs(eigs[2] - eigs[1]) < tol) -// { -// for(i = 0; i < 4; i++) ci[i] = 0.0; -// return; -// } - e1 = eigs[0]; e2 = eigs[1]; e3 = eigs[2]; @@ -438,12 +411,6 @@ static void fun_dos_case3(double en, double *eigs, double *ci) { double e1, e2, e3, e4; -// if(fabs(eigs[3] - eigs[2]) < tol) -// { -// for(i = 0; i < 4; i++) ci[i] = 0.0; -// return; -// } - e1 = eigs[0]; e2 = eigs[1]; e3 = eigs[2]; @@ -463,49 +430,53 @@ static void fun_dos_case3(double en, double *eigs, double *ci) static double F(double en, double e1, double e2, double e3, double e4) { - double complex s; + std::complex s; if(fabs(e1 - e3) > tol && fabs(e4 - e2) > tol) return (e1 - en) * (en - e2) / ((e1 - e3) * (e4 - e2)); else { +// Regularization to avoid division by zero s = fmin(fabs(e3 - e1), fabs(e4 - e2)); s /= 100.0; - s = fmax(s, 1.0e-20) * I; + s = fmax(std::abs(s), 1.0e-20) * I; - return creal((e1 - en + s) * (en - e2 + s) / ((e1 - e3 + s) * (e4 - e2 + s))); + return std::real((e1 - en + s) * (en - e2 + s) / ((e1 - e3 + s) * (e4 - e2 + s))); } } static double K2(double en, double e1, double e2, double e3) { - double complex s; + std::complex s; if(fabs(e1 - e3) > tol && fabs(e1 - e2) > tol) return (en - e1) / ((e2 - e1) * (e3 - e1)); else { +// Regularization to avoid division by zero s = fmin(fabs(e3 - e1), fabs(e1 - e2)); s /= 100.0; - s = fmax(s, 1.0e-20) * I; + s = fmax(std::abs(s), 1.0e-20) * I; - return creal((en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s))); + return std::real((en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s))); } } static double K1(double en, double e1, double e2) { - double complex s; + std::complex s; if(fabs(e1 - e2) > tol) return (e1 - en) / ((e2 - e1) * (e2 - e1)); else { +// Regularization to avoid division by zero s = fabs(e1 - e2); s /= 100.0; - s = fmax(s, 1.0e-20) * I; + s = fmax(std::abs(s), 1.0e-20) * I; - return creal((e1 - en + s) / ((e2 - e1 + s) * (e2 - e1 + s))); + return std::real((e1 - en + s) / ((e2 - e1 + s) * (e2 - e1 + s))); } } + diff --git a/c++/plovasp/atm/dos_tetra3d.h b/c++/plovasp/atm/dos_tetra3d.h deleted file mode 100644 index bc562e24..00000000 --- a/c++/plovasp/atm/dos_tetra3d.h +++ /dev/null @@ -1,42 +0,0 @@ -/******************************************************************************* -* -* This file is part of the ATM library. -* -* Copyright (C) 2010 by O. E. Peil -* -* TRIQS is free software: you can redistribute it and/or modify it under the -* terms of the GNU General Public License as published by the Free Software -* Foundation, either version 3 of the License, or (at your option) any later -* version. -* -* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY -* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more -* details. -* -* You should have received a copy of the GNU General Public License along with -* TRIQS. If not, see . -* -*******************************************************************************/ -#ifndef __C_DOS_TETRA3D_H__ -#define __C_DOS_TETRA3D_H__ - -#include -#include - -static PyObject *tetra_DOS3D(PyObject *self, PyObject *args); - -//void tet_dos3d(double en, double *eigk, int strd_eigk, -// npy_int64 *itt, int ntet, int *strd_itt, -// double *cti, int *strd_cti); -void tet_dos3d(double en, PyArrayObject *py_eigk, - PyArrayObject *py_itt, int ntet, - PyArrayObject *py_cti); -int dos_corner_weights(double en, double *eigs, int *inds, - double *ci); - -int dos_reorder(double en, double *e, int *inds); - -static const double small = 2.5e-2, tol = 1e-8; - -#endif diff --git a/c++/plovasp/atm/dos_tetra3d.hpp b/c++/plovasp/atm/dos_tetra3d.hpp new file mode 100644 index 00000000..8b07d1dd --- /dev/null +++ b/c++/plovasp/atm/dos_tetra3d.hpp @@ -0,0 +1,46 @@ +/******************************************************************************* +* +* This file is part of the ATM library. +* +* Copyright (C) 2010 by O. E. Peil +* +* TRIQS is free software: you can redistribute it and/or modify it under the +* terms of the GNU General Public License as published by the Free Software +* Foundation, either version 3 of the License, or (at your option) any later +* version. +* +* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY +* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +* details. +* +* You should have received a copy of the GNU General Public License along with +* TRIQS. If not, see . +* +*******************************************************************************/ +#pragma once + +#include + +using triqs::arrays::array; +using triqs::arrays::array_view; + +/// DOS of a band by analytical tetrahedron method +/// +/// Returns corner weights for all tetrahedra for a given band and real energy. +array_view +dos_tetra_weights_3d(array_view eigk, /// Band energies for each k-point + double en, /// Energy at which DOS weights are to be calculated + array_view itt /// Tetrahedra defined by k-point indices +); +//array +//dos_tetra_weights_3d(array eigk, /// Band energies for each k-point +// double e, /// Energy at which DOS weights are to be calculated +// array itt /// Tetrahedra defined by k-point indices +//); +int dos_corner_weights(double en, double *eigs, int *inds, double *ci); +int dos_tet_weights(double en, double *eigs, int *inds, double *ct); +int dos_reorder(double en, double *e, int *inds); + + +