diff --git a/c/vasp/atm/.gitignore b/c/vasp/atm/.gitignore new file mode 100644 index 00000000..675673c3 --- /dev/null +++ b/c/vasp/atm/.gitignore @@ -0,0 +1,5 @@ +makefile +Makefile +*.so +*.o +*.pyc diff --git a/c/vasp/atm/__init__.py b/c/vasp/atm/__init__.py new file mode 100644 index 00000000..139597f9 --- /dev/null +++ b/c/vasp/atm/__init__.py @@ -0,0 +1,2 @@ + + diff --git a/c/vasp/atm/argsort.c b/c/vasp/atm/argsort.c new file mode 100644 index 00000000..40347472 --- /dev/null +++ b/c/vasp/atm/argsort.c @@ -0,0 +1,47 @@ +#include + +int cmp(const void *a, const void *b) +{ + return (**(const double **)a) < (**(const double **)b) ? -1 : 1; +} + +int icmp(const void *a, const void *b) +{ + return (**(const int **)a) - (**(const int **)b); +} + +void argsort(double *arr, int *inds, double **ptrs, const int n) +{ + int i; + + for (i=0; i < n; i++) + { + ptrs[i] = arr + i; + } + + qsort(ptrs, n, sizeof(double *), cmp); + + for (i=0; i < n; i++) + { + inds[i] = (int)(ptrs[i] - arr); + } +} + +void iargsort(int *iarr, int *inds, int **ptrs, const int n) +{ + int i; + + for (i=0; i < n; i++) + { + ptrs[i] = iarr + i; + } + + qsort(ptrs, n, sizeof(int *), icmp); + + for (i=0; i < n; i++) + { + inds[i] = (int)(ptrs[i] - iarr); + } +} + + diff --git a/c/vasp/atm/argsort.h b/c/vasp/atm/argsort.h new file mode 100644 index 00000000..6bd581ca --- /dev/null +++ b/c/vasp/atm/argsort.h @@ -0,0 +1,7 @@ +#ifndef __ARGSORT_H__ +#define __ARGSORT_H__ + +void argsort(double *arr, int *inds, double **ptrs, const int n); +void iargsort(int *iarr, int *inds, int **ptrs, const int n); + +#endif diff --git a/c/vasp/atm/dos_tetra3d.c b/c/vasp/atm/dos_tetra3d.c new file mode 100644 index 00000000..baa4cdef --- /dev/null +++ b/c/vasp/atm/dos_tetra3d.c @@ -0,0 +1,464 @@ +#include + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION + +#include +#include +//#include "tetra.h" +#include "dos_tetra3d.h" + +/*************************************************** + + Analytical tetrahedron method as described in +Lambin et al., PRB 29, 6, 3430 (1984). + +***************************************************/ + +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); + +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); + +// +// 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 *eigk, *cti; // C-pointer to 'eigk' and 'rti' arrays + double e, en, eigs[4], ri[4]; + double et[4]; // Real part of eigenvalues + int strd_itt[2], strd_eigk, strd_cti[2], it_cti; // Strides + int inds[4], ntet; // Indices of sorted 'eigs'; number of tetrahedra +// Auxiliary variables and loop indices + int nd, it, ik, i, j, k, flag; + double *ptrs[4], e1e2, e2e3; + +// 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; + if (nd != 1) + { + PyErr_SetString(PyExc_ValueError, " Array 'eigk' must be 1D"); + return NULL; + } + + nd = py_itt->nd; + if (nd != 2) + { + PyErr_SetString(PyExc_ValueError, " Array 'itt' must be 2D"); + return NULL; + } + + nd = py_itt->dimensions[0]; + if (nd != 5) + { + PyErr_SetString(PyExc_ValueError, + " The first dimension of 'itt' must be equal to 5"); + return NULL; + } + + ntet = (int) py_itt->dimensions[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(4 * ntet * sizeof(double)); + dims = (npy_intp *)malloc(2 * sizeof(npy_intp)); + + dims[0] = 4; // Magic number: the number of tetrahedron corneres + 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 + +// +// Main function +// + tet_dos3d(e, eigk, strd_eigk, itt, ntet, strd_itt, cti, strd_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) +{ + double eigs[4], ci[4]; + int i, j, it, ik, it_cti, inds[4], flag; +// **** DEBUG + double ct, ci_sum; + +// 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]]; + eigs[i - 1] = eigk[ik * strd_eigk]; + } + +// corner weights for a single tetrahedron + dos_corner_weights(en, eigs, inds, ci); + for(i = 0, ci_sum = 0.0; i < 4; i++) + ci_sum += ci[i]; + + j = dos_tet_weights(en, eigs, inds, &ct); + if(fabs(ct - ci_sum) > tol) + { + printf(" *** Error in weights: it = %d, flag = %d, en = %lf", it, j, en); + for(i = 0; i < 4; i++) + printf(", e[%d] = %lf", i, eigs[i]); + printf(", c_diff = %le\n", fabs(ct - ci_sum)); + return; + } +// 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)); + +// 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++) + { + j = inds[i] * strd_cti[0] + it_cti; + cti[j] = ci[i]; + } + + } // it = 1, ntet +} + +int dos_corner_weights(double en, double *eigs, int *inds, + double *ci) +{ + int flag, i; + flag = dos_reorder(en, eigs, inds); + + switch(flag) + { +// E1 <= E <= E2 + case 1: + fun_dos_case1(en, eigs, ci); + break; + +// E2 <= E <= E3 + case 2: + fun_dos_case2(en, eigs, ci); + break; + +// E3 <= E <= E4 + case 3: + fun_dos_case3(en, eigs, ci); + break; + +// E < E1 || E4 < E + case 4: + case 5: + for(i = 0; i < 4; i++) ci[i] = 0.0; + break; + +// E1 == E4 == E + case 6: + for(i = 0; i < 4; i++) ci[i] = 0.25; + break; + } + + return flag; +} + +int dos_tet_weights(double en, double *eigs, int *inds, + double *ct) +{ + double e1, e2, e3, e4; + double complex s; + int flag, i; + flag = dos_reorder(en, eigs, inds); + + e1 = eigs[0]; + e2 = eigs[1]; + e3 = eigs[2]; + e4 = eigs[3]; + + switch(flag) + { +// E1 <= E <= E2 + case 1: + if(fabs(e2 - e1) > tol && fabs(e3 - e1) > tol && fabs(e4 - e1) > tol) + *ct = 3.0 * (en - e1) * (en - e1) / ((e2 - e1) * (e3 - e1) * (e4 - e1)); + else + { + 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; + + *ct = 3.0 * creal((en - e1 + s) * (en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s) * (e4 - e1 + s))); + } + + break; + +// E2 <= E <= E3 + case 2: + if(fabs(e4 - e2) > tol && fabs(e3 - e2) > tol && fabs(e4 - e1) > tol && fabs(e3 - e1) > tol) + *ct = 3.0 * ( + (e3 - en) * (en - e2) / ((e4 - e2) * (e3 - e2) * (e3 - e1)) + + (e4 - en) * (en - e1) / ((e4 - e1) * (e4 - e2) * (e3 - e1))); + else + { + s = fmin(fabs(e3 - e2), fabs(e3 - e1)); + s = fmin(fabs(s), fabs(e4 - e1)); + s = fmin(fabs(s), fabs(e4 - e2)); + s /= 100.0; + s = fmax(s, 1.0e-20) * I; + + *ct = 3.0 * creal(( + (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)))); + } + break; + +// E3 <= E <= E4 + case 3: + if(fabs(e4 - e2) > tol && fabs(e4 - e3) > tol && fabs(e4 - e1) > tol) + *ct = 3.0 * (e4 - en) * (e4 - en) / ((e4 - e1) * (e4 - e2) * (e4 - e3)); + else + { + 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; + + *ct = 3.0 * creal((e4 - en + s) * (e4 - en + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e4 - e3 + s))); + } + + break; + +// E < E1 || E4 < E + case 4: + case 5: + *ct = 0.0; + break; + +// E1 == E4 == E + case 6: + *ct = 1.0; + break; + } + + return flag; +} + +int dos_reorder(double en, double *e, int *inds) +{ + double *ptrs[4], e_tmp[4]; + int i; + + for(i = 0; i < 4; i++) + e_tmp[i] = e[i]; + + argsort(e_tmp, inds, ptrs, 4); + + for(i = 0; i < 4; i++) + e[i] = e_tmp[inds[i]]; + + if((e[0] <= en && en <= e[3]) && fabs(e[3] - e[0]) < tol) return 6; + if(e[0] <= en && en <= e[1]) return 1; + if(e[1] <= en && en <= e[2]) return 2; + if(e[2] <= en && en <= e[3]) return 3; + if(en < e[0]) return 4; + if(e[3] < en) return 5; + return -1; +} + +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]; + e3 = eigs[2]; + e4 = eigs[3]; + + ci[0] = K2(en, e1, e2, e4) * F(en, e2, e1, e1, e3) + + K2(en, e1, e2, e3) * F(en, e3, e1, e1, e4) + + K2(en, e1, e3, e4) * F(en, e4, e1, e1, e2); + + ci[1] = -K1(en, e1, e2) * F(en, e1, e1, e3, e4); + + ci[2] = -K1(en, e1, e3) * F(en, e1, e1, e2, e4); + + ci[3] = -K1(en, e1, e4) * F(en, e1, e1, e2, e3); +} + +static void fun_dos_case2(double en, double *eigs, double *ci) +{ + double e1, e2, e3, e4; + int i; + +// 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]; + e4 = eigs[3]; + + ci[0] = 0.5 * (K1(en, e3, e1) * ( + F(en, e3, e2, e2, e4) + + F(en, e4, e1, e2, e4) + + F(en, e3, e1, e2, e4)) + + K1(en, e4, e1) * ( + F(en, e4, e1, e2, e3) + + F(en, e4, e2, e2, e3) + + F(en, e3, e1, e2, e3))); + + ci[1] = 0.5 * (K1(en, e3, e2) * ( + F(en, e3, e2, e1, e4) + + F(en, e4, e2, e1, e4) + + F(en, e3, e1, e1, e4)) + + K1(en, e4, e2) * ( + F(en, e3, e2, e1, e3) + + F(en, e4, e1, e1, e3) + + F(en, e4, e2, e1, e3))); + + ci[2] = 0.5 * (-K1(en, e2, e3) * ( + F(en, e3, e2, e1, e4) + + F(en, e4, e2, e1, e4) + + F(en, e3, e1, e1, e4)) - + K1(en, e1, e3) * ( + F(en, e3, e2, e2, e4) + + F(en, e4, e1, e2, e4) + + F(en, e3, e1, e2, e4))); + + ci[3] = 0.5 * (-K1(en, e2, e4) * ( + F(en, e3, e2, e1, e3) + + F(en, e4, e1, e1, e3) + + F(en, e4, e2, e1, e3)) - + K1(en, e1, e4) * ( + F(en, e4, e1, e2, e3) + + F(en, e4, e2, e2, e3) + + F(en, e3, e1, e2, e3))); +} + +static void fun_dos_case3(double en, double *eigs, double *ci) +{ + double e1, e2, e3, e4; + int i; + +// 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]; + e4 = eigs[3]; + + ci[0] = K1(en, e4, e1) * F(en, e4, e4, e2, e3); + + ci[1] = K1(en, e4, e2) * F(en, e4, e4, e1, e3); + + ci[2] = K1(en, e4, e3) * F(en, e4, e4, e1, e2); + + ci[3] = -K2(en, e4, e3, e1) * F(en, e4, e3, e2, e4) - + K2(en, e4, e2, e3) * F(en, e4, e2, e1, e4) - + K2(en, e4, e1, e2) * F(en, e4, e1, e3, e4); + +} + +static double F(double en, double e1, double e2, double e3, double e4) +{ + double complex s; + + if(fabs(e1 - e3) > tol && fabs(e4 - e2) > tol) + return (e1 - en) * (en - e2) / ((e1 - e3) * (e4 - e2)); + else + { + s = fmin(fabs(e3 - e1), fabs(e4 - e2)); + s /= 100.0; + s = fmax(s, 1.0e-20) * I; + + return creal((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; + + if(fabs(e1 - e3) > tol && fabs(e1 - e2) > tol) + return (en - e1) / ((e2 - e1) * (e3 - e1)); + else + { + s = fmin(fabs(e3 - e1), fabs(e1 - e2)); + s /= 100.0; + s = fmax(s, 1.0e-20) * I; + + return creal((en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s))); + } +} + +static double K1(double en, double e1, double e2) +{ + double complex s; + + if(fabs(e1 - e2) > tol) + return (e1 - en) / ((e2 - e1) * (e2 - e1)); + else + { + s = fabs(e1 - e2); + s /= 100.0; + s = fmax(s, 1.0e-20) * I; + + return creal((e1 - en + s) / ((e2 - e1 + s) * (e2 - e1 + s))); + } +} + diff --git a/c/vasp/atm/dos_tetra3d.h b/c/vasp/atm/dos_tetra3d.h new file mode 100644 index 00000000..f0f714a4 --- /dev/null +++ b/c/vasp/atm/dos_tetra3d.h @@ -0,0 +1,18 @@ +#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); +int dos_corner_weights(double en, double *eigs, int *inds, + double *ci); + +int dos_reorder(double en, double *e, int *inds); + + +#endif diff --git a/c/vasp/atm/setup.py b/c/vasp/atm/setup.py new file mode 100644 index 00000000..49784320 --- /dev/null +++ b/c/vasp/atm/setup.py @@ -0,0 +1,9 @@ + +from distutils.core import setup, Extension +import numpy + +c_atm_dos_mod = Extension('c_atm_dos', sources=['dos_tetra3d.c', 'argsort.c'], + include_dirs=[numpy.get_include()]) + +setup(ext_modules=[c_atm_dos_mod]) +