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); + + +