[ase-users] Phonons example

[Horsley, Matthew A.]

Hi,
I am new to ase, and I am using the Phonon example (I copied it below). In the example, there is a parameter N, which controls the size of the supercell. In the example, N=7, and the code runs fine. However, if I change N to something larger or smaller (I tried N=5,6,8), the code crashes with an error message (see below). Is there a fix for this?
Thanks!
Matt

Error msg:
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-18-03ff49a02177> in <module>()
     14
     15 # Read forces and assemble the dynamical matrix
---> 16 ph.read(acoustic=True)
     17
     18 # High-symmetry points in the Brillouin zone

ase\phonons.pyc in read(self, method, symmetrize, acoustic, cutoff, born, **kwargs)
    404
    405                 # Slice out included atoms
--> 406                 C_Nav = C_av.reshape((N, len(self.atoms), 3))[:, self.indices]
    407                 index = 3 * i + j
    408                 C_xNav[index] = C_Nav

ValueError: cannot reshape array of size 1029 into shape (125,1,3)

Phonon example:

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.dft.kpoints import ibz_points, bandpath
from ase.phonons import Phonons

# Setup crystal and EMT calculator
atoms = bulk('Al', 'fcc', a=4.05)
calc = EMT()

# Phonon calculator
N = 7
ph = Phonons(atoms, calc, supercell=(N, N, N), delta=0.05)
ph.run()

# Read forces and assemble the dynamical matrix
ph.read(acoustic=True)

# High-symmetry points in the Brillouin zone
points = ibz_points['fcc']
G = points['Gamma']
X = points['X']
W = points['W']
K = points['K']
L = points['L']
U = points['U']

point_names = ['$\Gamma$', 'X', 'U', 'L', '$\Gamma$', 'K']
path = [G, X, U, L, G, K]

# Band structure in meV
path_kc, q, Q = bandpath(path, atoms.cell, 100)
omega_kn = 1000 * ph.band_structure(path_kc)

# Calculate phonon DOS
omega_e, dos_e = ph.dos(kpts=(50, 50, 50), npts=5000, delta=5e-4)
omega_e *= 1000

# Plot the band structure and DOS
import matplotlib.pyplot as plt
plt.figure(1, (8, 6))
plt.axes([.1, .07, .67, .85])
for n in range(len(omega_kn[0])):
    omega_n = omega_kn[:, n]
    plt.plot(q, omega_n, 'k-', lw=2)

plt.xticks(Q, point_names, fontsize=18)
plt.yticks(fontsize=18)
plt.xlim(q[0], q[-1])
plt.ylabel("Frequency ($\mathrm{meV}$)", fontsize=22)
plt.grid('on')

plt.axes([.8, .07, .17, .85])
plt.fill_between(dos_e, omega_e, y2=0, color='lightgrey', edgecolor='k', lw=1)
plt.ylim(0, 35)
plt.xticks([], [])
plt.yticks([], [])
plt.xlabel("DOS", fontsize=18)
plt.show()

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20180412/e1c6b9a8/attachment-0001.html>