Skip to content

Commit

Permalink
Add adaptive broadening (#131)
Browse files Browse the repository at this point in the history 
* Rename sc_offset for consistency with cell_origins

* Add initial mode gradients calculation

* Add initial per-mode broadening to DOS

* Fix adaptive dos indexing with small widths

* Always calculate Cart vectors

* Fix adaptive dos bin summation and bin indices

* Add partial dyn mat gradients calculation in C

* Add initial mode gradients calc in C (needs debug)

* Explicitly set dmat_grad to NULL

* Fix mode gradients calculation

- Use correct cartesian cell origins data type (use double rather than
  int)
- Use correct phase convention
- Correctly conj eigenvectors
- Fix weighting calculation
- Temporarily calculate for whole matrix (rather than just upper
  triangular)

* Fix mode grad units and adaptive dos slicing

* Add optional use of pdf in adaptive broadening

* Add function to read from phonon_dos

* Add tests for .phonon_dos file read

* Add tests for mode gradients calculation

* Enable read of DOS from phonon_dos and add tests

* Improve calculate_dos

 - Allow specification of either bin edges or centres with are_bin_edges
   (bin centres are preferable for adaptive broadening, whereas bin
   edges are preferable for binning with histogram)
 - Remove gaussian implementation (potential high memory useage and
   doesn't work as implemented)
 - Set minimum mode width

* Fix Python mode gradients calculation

* Remove duplicated 2 in mode widths scaling

* Don't calculate modg in Python unless required

* Optimise Python mode gradient calc

* Optimise Python dmat grad calc

- Make use of einsum
- Precalculate all combined (sc image + cell) origins

* Replace tabs with spacing in C files

* Ensure all_cell_origins is the correct type

* Update tests:

 - Allow test of pdos in test_castep_reader
 - Increase mode gradients tolerance for calculate_qpoint_frequencies
 - Remove mode-widths scale from adaptive DOS tests

* Use const bool in C to control calc of gradients

Hopefully this will allow better path prediction and reduce performance
loss

* Reduce number of loops calculating dyn mat in C

Only the upper triangular needs to be calculated, as both
the dynamical matrix and its gradients are Hermitian. When initially
implementing the gradients calculation, the whole matrix was calculated
for simplicity. Now that it has been tested, only the upper triangular
needs to be calculated, and the rest filled in. This brings the
dynamical matrix calculation performance back into line with what is
implemented in master.

* Output mode widths rather than gradients

* Update adaptive DOS test data

Results changed slightly in commits 9ff7335 and 5075a83 due to
improving the method of conversion from mode gradients to mode
widths, from using an arbitrary mode scaling to estimating the
q-spacing

* Update mode gradients -> mode widths test files

* Update tests to show DOS only changed by a scale factor

* Regenerate DOS data

* Regenerate DOS script data with new scaling

* Regenerate derived DOS test data

* Workaround for old Numpy vers reading phonon_dos

* Update mode widths docstring

* Add tests for bin edges/centres in calculate_dos

* Add adaptive DOS example to user docs

* Initial update to formulation in docstring based on feedback

* Add metadata to .phonon_dos reads

* Regenerate powder test data for new DOS implementation

It has been checked that the new results differ only by a scale
factor, as expected

* Remove are_bin_edges kwarg from calculate_dos

* Add mode gradients to widths utility function

This can be useful for converting mode gradients read from a CASTEP
.phonon_dos file to mode widths for use in calculate_dos

* Add adaptive broadening to euphonic-dos tool

* Fix mode_widths calc for low pint/Numpy versions

* Simplify conversion to mode gradients

* Small updates from review comments

- Add option to set minimum adaptive broadening width to calculate_dos
- Add comment to docstring of mode_gradients_to_widths about symmetry
  reduction giving inaccurate widths
- Remove commented out code

* Update from_castep_phonon_dos user docs

* Update CHANGELOG.rst
  • Loading branch information
@rebeccafair
rebeccafair authored Apr 15, 2021
1 parent 488af52 commit c060807
Show file tree
Hide file tree
Showing 48 changed files with 733,571 additions and 4,515 deletions.
12 changes: 12 additions & 0 deletions CHANGELOG.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,15 @@
``get_dispersion``
- ``Spectrum1D``, ``Spectrum1DCollection`` and ``Spectrum2D`` objects
have a new ``metadata`` attribute, see their docstrings for details
- Euphonic can now read DOS/PDOS from CASTEP .phonon_dos files with
``Spectrum1D.from_castep_phonon_dos`` and
``Spectrum1DCollection.from_castep_phonon_dos``
- **Adaptive broadening** is now available for DOS, which can obtain a
more representative DOS than standard fixed-width broadening. See
`the docs <https://euphonic.readthedocs.io/en/latest/dos.html#adaptive-broadening>`_
for details
- Adaptive broadening can be used in the ``euphonic-dos`` tool with the
``--adaptive`` argument

- Improvements:

Expand Down Expand Up @@ -62,6 +71,9 @@
- Default Monkhorst-Pack meshes (i.e. [6, 6, 6] in ``euphonic-dos``
and [20, 20, 20] in ``sample_sphere_structure_factor()``) have
been replaced by default grid-spacing values.

- The scaling of density of states has changed, due to a change
in implementation

`v0.4.0 <https://github.com/pace-neutrons/Euphonic/compare/v0.3.2...v0.4.0>`_
----------
Expand Down
38 changes: 33 additions & 5 deletions c/_euphonic.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define NPY_NO_DEPRECATED_API NPY_1_9_API_VERSION
#define PY_ARRAY_UNIQUE_SYMBOL EUPHONIC_NPY_ARRAY_API
#include <string.h>
#include <stdbool.h>
#include <omp.h>
#include <Python.h>
#include <numpy/arrayobject.h>
Expand All @@ -25,6 +26,8 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
PyArrayObject *py_dmat_weighting;
PyArrayObject *py_evals;
PyArrayObject *py_dmats;
PyArrayObject *py_modegs;
PyArrayObject *py_all_ogs_cart;
int dipole;
int reciprocal_asr;
int splitting;
Expand Down Expand Up @@ -62,6 +65,8 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
double *dmat_weighting;
double *evals;
double *dmats;
double *modegs;
double *all_ogs_cart;
int *n_sc_ims;
int *sc_im_idx;
int *cell_ogs;
Expand All @@ -79,6 +84,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
int n_cells;
int n_rqpts;
int dmats_len;
int modegs_len;
int n_split_qpts;
int q, i, qpos;
int max_ims;
Expand All @@ -88,7 +94,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
int n_gvecs;

// Parse inputs
if (!PyArg_ParseTuple(args, "OO!O!O!O!O!O!O!O!O!iiiO!O!is",
if (!PyArg_ParseTuple(args, "OO!O!O!O!O!O!O!O!O!iiiO!O!O!O!is",
&py_idata,
&PyArray_Type, &py_cell_vec,
&PyArray_Type, &py_recip_vec,
Expand All @@ -104,6 +110,8 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
&splitting,
&PyArray_Type, &py_evals,
&PyArray_Type, &py_dmats,
&PyArray_Type, &py_modegs,
&PyArray_Type, &py_all_ogs_cart,
&n_threads,
&scipy_dir)) {
return NULL;
Expand Down Expand Up @@ -152,13 +160,16 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
dmat_weighting = (double*) PyArray_DATA(py_dmat_weighting);
evals = (double*) PyArray_DATA(py_evals);
dmats = (double*) PyArray_DATA(py_dmats);
modegs = (double*) PyArray_DATA(py_modegs);
all_ogs_cart = (double*) PyArray_DATA(py_all_ogs_cart);
n_sc_ims = (int*) PyArray_DATA(py_n_sc_ims);
sc_im_idx = (int*) PyArray_DATA(py_sc_im_idx);
cell_ogs = (int*) PyArray_DATA(py_cell_ogs);
n_cells = PyArray_DIMS(py_fc)[0];
n_rqpts = PyArray_DIMS(py_rqpts)[0];
n_split_qpts = PyArray_DIMS(py_split_idx)[0];
dmats_len = PyArray_DIMS(py_dmats)[0];
modegs_len = PyArray_DIMS(py_modegs)[0];
max_ims = PyArray_DIMS(py_sc_im_idx)[3];
dmat_elems = 2*9*n_atoms*n_atoms;
if (dipole) {
Expand All @@ -185,7 +196,8 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
omp_set_num_threads(n_threads);
#pragma omp parallel
{
double *corr, *dmat_per_q;
const bool calc_dmat_grad = (modegs_len > 0) ? true : false;
double *corr, *dmat_per_q, *dmat_grad;
if (dipole) {
corr = (double*) malloc(dmat_elems*sizeof(double));
}
Expand All @@ -195,9 +207,16 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
if (dmats_len == 0) {
dmat_per_q = (double*) malloc(dmat_elems*sizeof(double));
}
// If space for the mode gradients has been allocated, assume they
// should be calculated and allocate space for dyn mat gradients
if (calc_dmat_grad) {
dmat_grad = (double*) malloc(3*dmat_elems*sizeof(double));
} else {
dmat_grad = NULL;
}
#pragma omp for
for (q = 0; q < n_rqpts; q++) {
double *qpt, *dmat, *eval;
double *qpt, *dmat, *eval, *modeg;
qpt = (rqpts + 3*q);
eval = (evals + q*3*n_atoms);

Expand All @@ -206,8 +225,12 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
} else {
dmat = (dmats + q*dmat_elems);
}
if (calc_dmat_grad) {
modeg = (modegs + q*3*n_atoms);
}
calculate_dyn_mat_at_q(qpt, n_atoms, n_cells, max_ims, n_sc_ims,
sc_im_idx, cell_ogs, sc_ogs, fc, dmat);
sc_im_idx, cell_ogs, sc_ogs, fc, all_ogs_cart, calc_dmat_grad,
dmat, dmat_grad);

if (dipole) {
calculate_dipole_correction(qpt, n_atoms, cell_vec, recip_vec,
Expand Down Expand Up @@ -236,9 +259,14 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
add_arrays(dmat_elems, corr, dmat);
}

mass_weight_dyn_mat(dmat_weighting, n_atoms, dmat);
mass_weight_dyn_mat(dmat_weighting, n_atoms, 2, dmat);
diagonalise_dyn_mat_zheevd(n_atoms, qpt, dmat, eval, zheevd);
evals_to_freqs(n_atoms, eval);

if (calc_dmat_grad) {
mass_weight_dyn_mat(dmat_weighting, n_atoms, 6, dmat_grad);
calculate_mode_gradients(n_atoms, eval, dmat, dmat_grad, modeg);
}
}
}

Expand Down
115 changes: 101 additions & 14 deletions c/dyn_mat.c
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#define PY_SSIZE_T_CLEAN
#define NPY_NO_DEPRECATED_API NPY_1_9_API_VERSION
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
Expand All @@ -12,12 +13,15 @@
void calculate_dyn_mat_at_q(const double *qpt, const int n_atoms,
const int n_cells, const int max_images, const int *n_sc_images,
const int *sc_image_i, const int *cell_origins, const int *sc_origins,
const double *fc_mat, double *dyn_mat) {
const double *fc_mat, const double *all_origins_cart, const bool calc_dmat_grad,
double *dyn_mat, double *dmat_grad) {

int i, j, n, nc, k, sc, ii, jj, idx;
int i, j, n, nc, k, sc, ii, jj, sc_img_idx, idx, idx_t;
double qdotr;
double phase_r;
double phase_i;
double rcart;
double phase[2];
double phase_sum[2];
double rcart_sum[6];

// Note: C calculated dynamical matrix uses e^-i(q.r) convention, whereas
// Python uses the e^i(q.r) convention. This differing convention is used
Expand All @@ -31,28 +35,69 @@ void calculate_dyn_mat_at_q(const double *qpt, const int n_atoms,
int s_fc = 9*n_atoms*n_atoms; // For fc_mat

memset(dyn_mat, 0, 2*9*n_atoms*n_atoms*sizeof(double));
if (calc_dmat_grad) {
memset(dmat_grad, 0, 3*2*9*n_atoms*n_atoms*sizeof(double));
}
for (i = 0; i < n_atoms; i++) {
for (j = i; j < n_atoms; j++) {
for (nc = 0; nc < n_cells; nc++){
phase_r = 0;
phase_i = 0;
memset(phase_sum, 0, 2*sizeof(double));
memset(rcart_sum, 0, 6*sizeof(double));
// Calculate and sum phases for all images
for (n = 0; n < n_sc_images[nc*s_n[0] + i*s_n[1] + j]; n++) {
qdotr = 0;
sc = sc_image_i[nc*s_i[0] + i*s_i[1] + j*s_i[2] + n];
sc_img_idx = nc*s_i[0] + i*s_i[1] + j*s_i[2] + n;
sc = sc_image_i[sc_img_idx];
for (k = 0; k < 3; k++){
qdotr += qpt[k]*(sc_origins[3*sc + k] + cell_origins[3*nc + k]);
}
phase_r += cos(2*PI*qdotr);
phase_i -= sin(2*PI*qdotr);
phase[0] = cos(2*PI*qdotr);
phase[1] = -sin(2*PI*qdotr);
phase_sum[0] += phase[0];
phase_sum[1] += phase[1];
if (calc_dmat_grad) {
for (k = 0; k < 3; k++){
// Note: use cos + isin phase as dyn mat gradients aren't passed
// to a Fortran lib so we need to use the e^i(q.r) convention
rcart = all_origins_cart[3*sc_img_idx + k];
rcart_sum[2*k] += phase[1]*rcart; //Multiply phase by i: swap re and im
rcart_sum[2*k + 1] += phase[0]*rcart;
}
}
}
for (ii = 0; ii < 3; ii++){
for (jj = 0; jj < 3; jj++){
idx = (3*i+ii)*3*n_atoms + 3*j + jj;
// Real part
dyn_mat[2*idx] += phase_r*fc_mat[nc*s_fc + idx];
dyn_mat[2*idx] += phase_sum[0]*fc_mat[nc*s_fc + idx];
// Imaginary part
dyn_mat[2*idx + 1] += phase_i*fc_mat[nc*s_fc + idx];
dyn_mat[2*idx + 1] += phase_sum[1]*fc_mat[nc*s_fc + idx];
if (calc_dmat_grad) {
for (k = 0; k < 3; k++) {
// Real
dmat_grad[6*idx + 2*k] += rcart_sum[2*k]*fc_mat[nc*s_fc + idx];
// Imaginary
dmat_grad[6*idx + 2*k + 1] += rcart_sum[2*k + 1]*fc_mat[nc*s_fc + idx];
}
}
}
}
}
}
}

// Fill in lower triangular of dmat gradients - is Hermitian
if (calc_dmat_grad) {
for (i = 1; i < n_atoms; i++) {
for (j = 0; j < i; j++) {
for (ii = 0; ii < 3; ii++){
for (jj = 0; jj < 3; jj++){
idx = 6*((3*i+ii)*3*n_atoms + 3*j + jj);
idx_t = 6*((3*j+jj)*3*n_atoms + 3*i + ii);
for (k = 0; k < 3; k++) {
dmat_grad[idx + 2*k] = dmat_grad[idx_t + 2*k];
dmat_grad[idx + 2*k + 1] = -dmat_grad[idx_t + 2*k + 1];
}
}
}
}
Expand Down Expand Up @@ -280,12 +325,14 @@ void calculate_gamma_correction(const double q_dir[3], const int n_atoms,
}

void mass_weight_dyn_mat(const double* dyn_mat_weighting, const int n_atoms,
double* dyn_mat) {
const int repeats, double* dyn_mat) {

int i, j;
for (i = 0; i < 9*n_atoms*n_atoms; i++) {
for (j = 0; j < 2; j++) {
dyn_mat[2*i + j] *= dyn_mat_weighting[i];
for (j = 0; j < repeats; j++) {
// Repeats: how many elements of dyn_mat per dyn_mat_weighting
// As dyn_mat = complex and weighting = real, this is usually 2
dyn_mat[repeats*i + j] *= dyn_mat_weighting[i];
}
}
}
Expand Down Expand Up @@ -347,3 +394,43 @@ void evals_to_freqs(const int n_atoms, double *eigenvalues) {
eigenvalues[i] = tmp;
}
}

void calculate_mode_gradients(const int n_atoms, const double *evals,
const double *evecs, const double *dmat_grad, double *modeg) {
int n, i, j, a, b, k, grad_idx;
int n_modes = 3*n_atoms;
double grad_dot;
double grad_tmp[6];
double evec_mult_tmp[2];
double conj_tmp[2];
int mode_s = 2*3*n_atoms; //Eigenvector array stride

for (n = 0; n < n_modes; n++) {
memset(grad_tmp, 0, 6*sizeof(double));
for (i = 0; i < n_atoms; i++) {
for (a = 0; a < 3; a++) {
for (j = 0; j < n_atoms; j++) {
for (b = 0; b < 3; b++) {

for (k = 0; k < 3; k++) {
cmult_conj((evecs + (n*mode_s + 6*j + 2*b)),
(evecs + (n*mode_s + 6*i + 2*a)),
evec_mult_tmp);
grad_idx = 3*(3*i + a)*mode_s + 3*(6*j + 2*b) + 2*k;
cmult((dmat_grad + grad_idx), evec_mult_tmp, conj_tmp);
grad_tmp[2*k] += conj_tmp[0];
grad_tmp[2*k + 1] += conj_tmp[1];
}
}
}
}
}
grad_dot = 0;
for (k = 0; k < 3; k++) {
cmult_conj((grad_tmp + 2*k), (grad_tmp + 2*k), conj_tmp);
grad_dot += conj_tmp[0];
}
modeg[n] = 0.5*sqrt(grad_dot)/evals[n];
}
}

8 changes: 6 additions & 2 deletions c/dyn_mat.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,8 @@
void calculate_dyn_mat_at_q(const double *qpt, const int n_atoms,
const int n_cells, const int max_ims, const int *n_sc_images,
const int *sc_image_i, const int *cell_origins, const int *sc_origins,
const double *fc_mat, double *dyn_mat);
const double *fc_mat, const double *all_origins_cart,
const bool calc_dmat_grad, double *dyn_mat, double *dmat_grad);

void calculate_dipole_correction(const double *qpt, const int n_atoms,
const double *cell_vec, const double *recip, const double *atom_r,
Expand All @@ -18,11 +19,14 @@ void calculate_gamma_correction(const double q_dir[3], const int n_atoms,
const double *dielectric, double *corr);

void mass_weight_dyn_mat(const double *dyn_mat_weighting, const int n_atoms,
double *dyn_mat);
const int repeats, double *dyn_mat);

int diagonalise_dyn_mat_zheevd(const int n_atoms, const double qpt[3],
double *dyn_mat, double *eigenvalues, ZheevdFunc zheevdptr);

void evals_to_freqs(const int n_atoms, double *eigenvalues);

void calculate_mode_gradients(const int n_atoms, const double *evals,
const double *evecs, const double *dmat_grad, double *modeg);

#endif
Loading

0 comments on commit c060807

Please sign in to comment.