[gpaw-users] Phonon band structure

[Atz Togo]

Dear GPAW-users,

I am a post-doc in Kyoto university in Japan. I work for phonopy.
The Si band structure may be calculated as follows:

from gpaw import GPAW
from ase import *
from ase.calculators import numeric_force
from phonopy import Phonopy
import numpy as np

# GPAW setting
a = 5.404
bulk = Atoms(symbols='Si8',
             positions=[(0, 0, 0),
                        (0, 0.5, 0.5),
                        (0.5, 0, 0.5),
                        (0.5, 0.5, 0),
                        (0.25, 0.25, 0.25),
                        (0.25, 0.75, 0.75),
                        (0.75, 0.25, 0.75),
                        (0.75, 0.75, 0.25)],
             pbc=True)
bulk.set_cell((a, a, a), scale_atoms=True)

n = 32
calc = GPAW(gpts=(n, n, n),
            nbands=8*3,
            width=0.01,
            kpts=(2, 2, 2),
            convergence={'eigenstates': 1e-7}
            )

# Phonopy pre-process
print "------"
print "Phonon"
print "------"
# 1st arg. is the input unit cell.
# 2nd arg. is the supercell lattice relative to the unit cell.
# 'distance' is the distance of displacements.
# Default symmetry tolerance is 1e-5 in fractional coordinates.
phonon = Phonopy(bulk, [[1,0,0],[0,1,0],[0,0,1]], distance=0.01)
symmetry = phonon.get_symmetry()
print "Space group:", symmetry.get_international_table()


phonon.print_displacements()
supercells = phonon.get_supercells_with_displacements()

# Force calculations
set_of_forces = []
for scell in supercells:
    cell = Atoms( symbols=scell.get_chemical_symbols(),
                  scaled_positions=scell.get_scaled_positions(),
                  cell=scell.get_cell(),
                  pbc=True )
    cell.set_calculator(calc)
    forces = cell.get_forces()
    drift_force = forces.sum(axis=0)
    print "        ---------------------------------"
    print "     ", "%11.5f"*3 % tuple(drift_force)
    # Simple translational invariance
    for force in forces:
        force -= drift_force / forces.shape[0]
    set_of_forces.append(forces)

# Phonopy post-process
# 1st arg. is a relative lattice to the input unit cell.
# 2nd arg. is bunch of the calculated forces.
phonon.set_post_process([[0,0.5,0.5],[0.5,0,0.5],[0.5,0.5,0]], set_of_forces)

# Phonon dispersion relation
bands = []

q_start  = np.array([0.5,0.5,0.0])
q_end    = np.array([0.0,0.0,0.0])
band = []
for i in range(51):
    band.append( q_start + ( q_end - q_start ) / 50 * i )

bands.append(band)

q_start  = np.array([0.0,0.0,0.0])
q_end    = np.array([0.5,0.0,0.0])
band = []
for i in range(51):
    band.append( q_start + ( q_end - q_start ) / 50 * i )

bands.append(band)

#*********************
# Matplotlib required
#*********************
print "\nPhonon dispersion:"
phonon.set_band_structure( bands )
phonon.plot_band_structure().show()





Togo

On Fri, Dec 3, 2010 at 7:02 PM, Kristen Kaasbjerg <kkaa at fysik.dtu.dk> wrote:
> If someone finds it useful then why not ?
> I should stress though, that the ASE implementation has features for
> calculating phonon dispersion and mode visualization with ag.
>
> By the way, frophon code seem to have been replaced by the Phononpy code:
> http://phonopy.sourceforge.net/
>
> Kristen
>
>> Hi,
>>
>> by looking at the latest abinit release one can see that they
>> "Print out VASP-style POSCAR and FORCES files, for use with PHON or
>> frophon codes for frozen phonon calculations."
>> http://www.abinit.org/documentation/helpfiles/for-v6.4/release_notes/release_notes.html
>>
>> Could this way be also interesting for us?
>> website of phon: http://www.homepages.ucl.ac.uk/~ucfbdxa/phon/
>>
>> Marcin
>>
>
> _______________________________________________
> gpaw-users mailing list
> gpaw-users at listserv.fysik.dtu.dk
> https://listserv.fysik.dtu.dk/mailman/listinfo/gpaw-users
>



-- 
Atsushi Togo
http://atztogo.users.sourceforge.net/
atz.togo at gmail.com