In reply to @tomasusi:matrix.org
shouldn't it be a convergence parameter?

In reply to @tomasusi:matrix.org
the reason is that the multislice simulation doesn't know anything about energy loss, so we need to have the "correct phonon energy" in the structures

In reply to @tomasusi:matrix.org
it's a clever way of avoiding describing quite theoretically demanding inelastic momentum-resolved scattering

In reply to @flokno:matrix.org
I'd actually be interested to learn more about this

In reply to @tomasusi:matrix.org
That's the original paper:
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.104.104301

In reply to @chronum:matrix.org
I'm fairly certain N corresponds to the supercell counting index.

In reply to @chronum:matrix.org
I don't, unfortunately. At that point you're likely better off looking at ASE's source code.

Can I get a little help with this? I attached my code and output, and can provide anything else, I can also make this a gitlab issue if that would be more appropriate.

I set up three carbon atoms ~1.5 angstroms away from each other in a line on the x axis, and I expect the dynamical matrix to be of the form of a 9x9 matrix: [[2F, -F, 0], [-F, 2F, F], [0, -F, 2F]] where F=[[f, 0, 0], [0,0,0], [0,0,0]] (image attached). However, using the Phonons class, the dynamical matrix gives me only F--instead of 2F--in the upper left and lower right.

The way I set up the phonon calculation is to use a 3x1x1 supercell of non-periodic unit cells. Could this have somehow caused the deviation from the theoretical calculation?

Also there were some non-zero diagonal elements on F, which is strange since all the atoms were perfectly positioned along the x-axis (although they were usually only about a tenth of the value of f). Is this expected numerical error even though the setup was perfectly x-aligned?

Code:

from ase.phonons import Phonons
from ase import Atoms
from ase.optimize import BFGS

from gpaw import GPAW, PW

import numpy as np

# Defice atom positions and cell
atoms = Atoms("C3", positions=[(1.7432107, 0.75, 0.75), (3, 0.75, .75), (4.2567893, 0.75, 0.75)], pbc=[0,0,0])
atoms.set_cell((4.5, 3, 3))
atoms.center()

# Set the calculator
calc = GPAW(mode=PW(340), xc="PBE", txt="C3.txt", symmetry="off", convergence={'density': 1e-6, 'forces':1e-4})
atoms.calc = calc
print(atoms.get_potential_energy())
dyn = BFGS(atoms)
dyn.run()
print(atoms.get_potential_energy())

# Phonon calculations
ph = Phonons(atoms, calc, supercell=(3, 1, 1))
ph.run()
# ph.read(acoustic=True)
ph.read()
ph.clean()

np.set_printoptions(precision=1)
print(ph.D_N[0])

Output:

Pot_eng before BFGS: -13.911701139494758
      Step     Time          Energy          fmax
BFGS:    0 22:54:07      -13.911701        6.614344
BFGS:    1 22:54:10      -14.066930        5.722778
BFGS:    2 22:54:14      -14.284913        0.563426
BFGS:    3 22:54:16      -14.287304        0.044011
Pot_eng after BFGS: -14.287304483684558
Positions after BFGS:
[[1.04779882 1.5        1.5       ]
 [2.25       1.5        1.5       ]
 [3.45220118 1.5        1.5       ]]
-----
ph.D_N[0]:
[[ 9.7  0.   0.  -8.3  0.   0.   0.   0.   0. ]
 [ 0.  -1.3  0.   0.   1.   0.   0.  -0.4  0. ]
 [ 0.   0.  -1.3  0.   0.   1.   0.   0.  -0.4]
 [-8.3  0.   0.  16.9  0.   0.  -8.3  0.   0. ]
 [ 0.   1.   0.   0.  -1.5  0.   0.   1.   0. ]
 [ 0.   0.   1.   0.   0.  -1.5  0.   0.   1. ]
 [ 0.   0.   0.  -8.3  0.   0.   9.7  0.   0. ]
 [ 0.  -0.4  0.   0.   1.   0.   0.  -1.3  0. ]
 [ 0.   0.  -0.4  0.   0.   1.   0.   0.  -1.3]]