Fix units and force constant expansion in Phonons.load_phonopy#121
Fix units and force constant expansion in Phonons.load_phonopy#121krystophny wants to merge 2 commits into
Conversation
load_phonopy assumed Bohr coordinates and Ry/au^2 force constants regardless of the calculator that produced the phonopy files, and expanded the compact FORCE_CONSTANTS format with a wrong atom mapping, so the spectrum was wrong away from Gamma. The function was disabled with a NotImplementedError and its test was skipped. Read the physical_unit block of phonopy.yaml and convert lengths and force constants to the internal Angstrom and Ry/bohr^2 conventions, covering the unit sets phonopy writes for its calculators. Expand compact force constants through the lattice translation that connects each supercell atom to its primitive representative (the reduced_to field). Re-enable the TestPhonopyInput test and check structure and frequencies against values computed independently with phonopy, for both the au/Ry fixture and a converted Angstrom/eV fixture. Addresses SSCHAcode/python-sscha#167
mesonepigreco
left a comment
There was a problem hiding this comment.
Thank you very much for the contribution.
The previous version of the phonopy loader was not working, and we realized this much later due to a different convention between phonopy and cellconstructor in ordering atoms. The problem is that if you exchange blocks of atoms while reconstructing the supercell, the eigenvalues of the dynamical matrix are unaffected, but the eigenvectors will be. This problem does not show up when just diagonalizing and checking the frequencies of phonons, but becomes evident when you try to run an SSCHA calculation or when you try to Fourier-interpolate the dynamical matrix (e.g., to plot the dispersion).
The current test seems to benchmark only the frequencies (eigenvalues); therefore, a bug like the previous one is invisible to it. Probably the best strategy would be to test by computing a dynamical matrix using phonopy on a simple system (it should have at least 2 atoms in the primitive cell to spot wrong atomic assignment, possibly the same type, and at least some q points, e.g. a 3x3x3 supercell to include some non perfectly symmetrical) and compare with the result obtained from Phonons.compute_phonons_finite_displacements, which we know has the correct atomic assignment.
The existing TestPhonopyInput fixture has one atom per primitive cell and checks only the supercell frequencies, so it cannot detect a wrong atom assignment: permuting atoms leaves the eigenvalues unchanged and corrupts only the eigenvectors. Add TestPhonopyEigenvectors: two inequivalent same-species (Cu) atoms in a tetragonal cell, 3x3x3 supercell, compact FORCE_CONSTANTS so the reduced_to expansion is exercised. The test compares the full dynamical matrices against compute_phonons_finite_displacements, whose atom assignment is correct by construction. The reference is driven by an exact harmonic force law built from the same force constants and is committed precomputed, so the test needs only numpy; generate_reference.py reproduces it with phonopy/ase. The current loader matches the reference to 1.5e-6 Ry/bohr^2, while swapping the two atoms shifts the matrix by 9e-2 with identical frequencies.
|
Thanks, you're right: the Nb fixture has one atom per cell and the check was on eigenvalues, so it was blind to atom assignment. I added The current loader matches that reference to 1.5e-6 Ry/bohr^2. Swapping the two atoms leaves every frequency unchanged but shifts the matrix by 9e-2, so the test fails on exactly the assignment bug you described. The reference matrices are committed precomputed, so the test needs only numpy; |
Phonons.load_phonopy was disabled with a NotImplementedError because it returned wrong results. Two bugs:
Changes:
Addresses SSCHAcode/python-sscha#167
Verification
Before, on master, the new test fails (load_phonopy raises):
With the NotImplementedError bypassed, the old code gives a wrong spectrum for the bcc Nb 2x2x2 fixture. Reference computed independently with phonopy 2.x from the same files (cm^-1): 0 x3, 91.427 x12, 138.125 x6, 168.055 x6, 172.073 x6, 176.968 x6, 195.129 x6, 217.426 x3. Old code:
After the fix:
Full suite: