Fix phonopy dipole read (#241)

* Convert dielectric units

* Subtract supercell dipole corr from Phonopy fc

* Update NaCl fcs and add phonopy quartz test

* Update NAC tests/logic

- Don't require nac factor is in phonopy.yaml, even if born is present
- Don't require nac factor is in BORN file
- Raise KeyError if nac factor is in neither phonopy.yaml or BORN
- Only try to dipole correct force constants if born is present
- Add tests for the above combinations of nac in phonopy/BORN/neither
- Update NaCl_prim/phonopy.yaml to have nac_factor (otherwise causes
  error)

* Fix convert_fc_phases test

convert_fc_phases only reindexes the force constants matrix, and
doesn't remove the dipole from the force constants matrix (whereas
the dipole is now removed when reading from phonopy files). So we
cannot use the same NaCl_force_constants.json file as we use in the
ForceConstants.from_phonopy tests. So instead use a new file
NaCl_dipole_force_constants.json, which is reindexed, but doesn't
have the dipole removed (this is actually the same as the original
NaCl_force_constants.json file before the dipole fix was introduced).

* Update NaCl phonon modes test data

* Use NaCl-QE test case instead of quartz

* Add subtract_fc_dipole func to util.py

* Update NaCl script tests

* Remove erroneously added docstring

* Fix tests for old Pint versions

Older versions of Pint strip the Quantity with some Numpy operations
so use magnitudes and reapply units later to avoid this

* Small Codacy updates

* Convert calculate_dipole_correction to static method

* Read from dipole dict in C

* Cast force constants to real

* Move dipole func to ForceConstants

* Add test for from_long_ranged_dipole_fc

* Add user documentation

* Add docstring for static functions

* Rename NaCl_dipole to NaCl_long_ranged for clarity

* Changes from review

- Change user docs and docstring to be more clear about 'total' vs.
  'short-ranged' vs 'dipole-dipole' force constants
- Update instead of merge Born dict in from_phonopy

* Reorder force constants check

It's better to check easy-to-diagnose attributes first (e.g.
cell_origins, born charges) before testing the force constants
themselves, as they are much more manageable and will reveal
an issue much more easily than the effect this will have on the
force constants, which will just be small differences in a large
array of floating point numbers.

* Rename from_long_ranged_dipole_fc to from_total_fc_with_dipole

* Update changelog

CHANGELOG.rst

@@ -1,9 +1,18 @@
 `Unreleased <https://github.com/pace-neutrons/Euphonic/compare/v1.0.0...HEAD>`_
 ----------
 
+- New features:
+
+  - There is a new function ``ForceConstants.from_total_fc_with_dipole`` to allow
+    reading force constants from other programs which contain long-ranged
+    dipole-dipole interactions.
+
 - Bug fixes:
 
-  - Avoid occasional segmentation faults when using OpenBLAS, workaround for #191
+  - Avoid occasional segmentation faults when using OpenBLAS, workaround for
+    `#191 <https://github.com/pace-neutrons/Euphonic/issues/191>`_
+  - Correctly read force constants from Phonopy with dipole-dipole
+    interactions, see `#239 <https://github.com/pace-neutrons/Euphonic/issues/239>`_.
 
 `v1.0.0 <https://github.com/pace-neutrons/Euphonic/compare/v0.6.5...v1.0.0>`_
 ----------

c/_euphonic.c

@@ -40,6 +40,7 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
     // Extra vars only required if dipole = True
     PyArrayObject *py_born;
     PyArrayObject *py_dielectric;
+    PyDictObject *py_dipole_init_data;
     double lambda;
     PyArrayObject *py_H_ab;
     PyArrayObject *py_dipole_cells;
@@ -113,7 +114,6 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
                           &n_threads)) {
         return NULL;
     }
-
     // Load library functions
     // Load before calling PyObject_GetAttrString so we don't have to
     // py_DECREF if library loading fails
@@ -134,14 +134,19 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
             return NULL;
     }
     if (dipole) {
+        if (attr_from_pyobj(py_idata, "_dipole_init_data", &py_dipole_init_data)) {
+                PyErr_Format(PyExc_RuntimeError,
+                             "Failed to read dipole init data from object\n");
+                return NULL;
+        }
         if (attr_from_pyobj(py_idata, "_born", &py_born) ||
             attr_from_pyobj(py_idata, "_dielectric", &py_dielectric) ||
-            double_from_pyobj(py_idata, "_lambda", &lambda) ||
-            attr_from_pyobj(py_idata, "_H_ab", &py_H_ab) ||
-            attr_from_pyobj(py_idata, "_cells", &py_dipole_cells) ||
-            attr_from_pyobj(py_idata, "_gvec_phases", &py_gvec_phases) ||
-            attr_from_pyobj(py_idata, "_gvecs_cart", &py_gvecs_cart) ||
-            attr_from_pyobj(py_idata, "_dipole_q0", &py_dipole_q0)) {
+            double_from_pydict(py_dipole_init_data, "lambda", &lambda) ||
+            val_from_pydict(py_dipole_init_data, "H_ab", &py_H_ab) ||
+            val_from_pydict(py_dipole_init_data, "cells", &py_dipole_cells) ||
+            val_from_pydict(py_dipole_init_data, "gvec_phases", &py_gvec_phases) ||
+            val_from_pydict(py_dipole_init_data, "gvecs_cart", &py_gvecs_cart) ||
+            val_from_pydict(py_dipole_init_data, "dipole_q0", &py_dipole_q0)) {
                 PyErr_Format(PyExc_RuntimeError,
                              "Failed to read dipole attributes from object\n");
                 return NULL;
@@ -288,13 +293,11 @@ static PyObject *calculate_phonons(PyObject *self, PyObject *args) {
     Py_DECREF(py_cell_ogs);
     Py_DECREF(py_atom_r);
     if (dipole){
+        // Note PyDict_GetItemString returns "borrowed" ref so don't need
+        // to decref lambda, H_ab etc.
         Py_DECREF(py_born);
         Py_DECREF(py_dielectric);
-        Py_DECREF(py_H_ab);
-        Py_DECREF(py_dipole_cells);
-        Py_DECREF(py_gvec_phases);
-        Py_DECREF(py_gvecs_cart);
-        Py_DECREF(py_dipole_q0);
+        Py_DECREF(py_dipole_init_data);
     }
 
     return Py_None;

c/py_util.c

@@ -4,6 +4,33 @@
 #include <Python.h>
 #include <numpy/arrayobject.h>
 
+int val_from_pydict(PyDictObject *dict, const char *key, PyObject **result) {
+/* Given a PyDictObject and key, retrieve get the address of the value associated
+ * with that key, and alter the pointer pointed to by result to point to that
+ * address */
+    PyObject *tmp = PyDict_GetItemString(dict, key);
+    if (!tmp) {
+        printf("Couldn't retrieve %s from dict\n", key);
+        return 1;
+    }
+    *result = tmp;
+    return 0;
+}
+
+int double_from_pydict(PyDictObject *dict, const char *key, double *result) {
+/* Given a PyDictObject and key, read a double from that dict and store it in
+   the address pointed to by result */
+    PyObject *tmp;
+    val_from_pydict(dict, key, &tmp);
+    if (PyFloat_Check(tmp)) {
+        *result = (double) PyFloat_AsDouble(tmp);
+    } else {
+        printf("Incorrect type for %s\n", key);
+        return 1;
+    }
+    return 0;
+}
+
 int attr_from_pyobj(PyObject *obj, const char *attr_name, PyObject **result) {
 /* Given a PyObject and the name of one of its attributes, get the address of
  * that attribute, and alter the pointer pointed to by result to point to that

c/py_util.h

@@ -1,6 +1,8 @@
 #ifndef __py_util_H__
 #define __py_util_H__
 
+int val_from_pydict(PyDictObject *dict, const char *key, PyObject **result);
+int double_from_pydict(PyDictObject *dict, const char *key, double *result);
 int attr_from_pyobj(PyObject *obj, const char *attr_name, PyObject **result);
 int int_from_pyobj(PyObject *obj, const char *attr_name, int *result);
 int double_from_pyobj(PyObject *obj, const char *attr_name, double *result);

doc/source/force-constants.rst

@@ -87,6 +87,28 @@ This is a simple example, but the same applies to more realistic 3D structures.
 This transformation is done automatically when reading from Phonopy.
 If reading force constants from another program, there is a function to help with this transformation, see :ref:`Reading From Other Programs<fc_read_other_programs>`.
 
+**Dipole-dipole Force Constants**
+
+For materials without dipole-dipole interactions (non-zero Born effective charges) force constants typically decay as :math:`1/r^{5-7}`, but if dipole-dipole interactions are present, they only decay at :math:`1/r^3`.
+It is often not computationally practical to compute the force constants at such a distance.
+Fortunately there is an analytical expression for these dipole interactions, which can be applied to the dynamical matrix to produce the correct phonon frequencies and displacements.
+However, to avoid double counting some of the dipole-dipole interactions, the force constants matrix from which the dynamical matrix is calculated must not already include the short-ranged part of the dipole interactions.
+This is represented in equation 78 in `Gonze & Lee, 1997 <https://doi.org/10.1103/PhysRevB.55.10355>`_:
+
+.. math::
+
+  C = C^{SR} + C^{DD}
+
+Where :math:`C` is the 'total' force constants matrix containing both short-ranged interactions and long-ranged dipole interactions,
+:math:`C^{SR}` is the force constants matrix containing only the short-ranged interactions (i.e. no dipole interactions),
+and :math:`C^{DD}` is the force constants matrix containing only the long-ranged dipole interactions,
+and all matrices have the same shape and indexing and contain interactions to the same cut-off distance.
+To correctly apply the correction to the resulting dynamical matrix, :math:`C^{SR}` must be used, which is what Euphonic requires.
+Some codes (e.g. Phonopy) output :math:`C` so must be converted.
+This conversion is done automatically when reading from Phonopy, but if reading force constants from another program, there is a class method
+:py:meth:`ForceConstants.from_total_fc_with_dipole <euphonic.force_constants.ForceConstants.from_total_fc_with_dipole>`
+which can do this transformation. An example of this is shown in :ref:`Reading From Other Programs<fc_read_other_programs>`
+
 
 Reading From CASTEP
 -------------------
@@ -175,6 +197,8 @@ An example is shown below, assuming that the inputs are being loaded from existi
   # Create a ForceConstants object, using the Crystal object
   fc = ForceConstants(crystal, force_constants, sc_matrix, cell_origins)
 
+**Changing Phase Convention**
+
 If, as described in the :ref:`Force Constants Format<fc_format>` section, the source program uses the atomic coordinate phase convention, there may be some re-indexing required to get the force constants in the correct shape and form.
 There is a helper function :py:meth:`euphonic.util.convert_fc_phases <euphonic.util.convert_fc_phases>`
 which will convert a force constants of shape ``(n, N*n, 3, 3)``, to the shape required by Euphonic ``(N, 3*n, 3*n)``,
@@ -217,8 +241,35 @@ For more details see the function docstring. An example is below:
   # Create a ForceConstants object, using the Crystal object
   fc = ForceConstants(crystal, force_constants, sc_matrix, cell_origins)
 
-Once the force constants object has been created, it can be saved as a single portable JSON file using
-:py:meth:`ForceConstants.to_json_file <euphonic.force_constants.ForceConstants.to_json_file>`
+**Dipole-dipole Force Constants - Converting from Total to Short-ranged**
+
+If, as described in the :ref:`Force Constants Format<fc_format>` section, the force constants matrix contains the dipole interactions, these must be subtracted to produce the short-ranged force constants matrix before being used in Euphonic.
+This can be done using the class method
+:py:meth:`ForceConstants.from_total_fc_with_dipole <euphonic.force_constants.ForceConstants.from_total_fc_with_dipole>`.
+Note that the force constants must have a Euphonic-like shape before using this method. An example is below:
+
+.. code-block:: py
+
+  import numpy as np
+  from euphonic import ureg, Crystal, ForceConstants
+  from euphonic.util import convert_fc_phases
+
+  # Load arrays from files and apply any required units
+  cell_vectors = np.load('cell_vectors.npy')*ureg('angstrom')
+  atom_r = np.load('atom_r.npy')
+  atom_type = np.load('atom_type.npy')
+  atom_mass = np.load('atom_mass.npy')*ureg('amu')
+  force_constants = np.load('force_constants.npy')*ureg('meV/(angstrom**2)')
+  sc_matrix = np.load('sc_matrix.npy')
+  cell_origins = np.load('cell_origins.npy')
+  born = np.load('born.npy')*ureg('e')
+  dielectric = np.load('dielectric.npy')*ureg('e**2/(bohr*hartree)')
+
+  # Create a Crystal object
+  crystal = Crystal(cell_vectors, atom_r, atom_type, atom_mass)
+  # Create a ForceConstants object from the long-ranged force constants
+  fc = ForceConstants.from_total_fc_with_dipole(
+    crystal, force_constants, sc_matrix, cell_origins, born, dielectric)
 
 
 Calculating Phonon Frequencies and Eigenvectors
@@ -265,4 +316,4 @@ Docstring
 
 .. autoclass:: euphonic.force_constants.ForceConstants
    :members:
-   :exclude-members: force_constants, born, dielectric
+   :exclude-members: force_constants, born, dielectric