Add files via upload

TP-Vibrations/Al_ase_phonons.py

@@ -0,0 +1,48 @@
+from ase.build import bulk
+from ase.calculators.emt import EMT
+from ase.phonons import Phonons
+from ipywidgets import FloatSlider, IntSlider, interact, interact_manual, fixed
+
+def calculate_phonons(x):
+    # Setup crystal and EMT calculator
+    atoms = bulk('Al', 'fcc', a=x)#4.05)
+
+    # Phonon calculator
+    N = 7
+    ph = Phonons(atoms, EMT(), supercell=(N, N, N), delta=0.05)
+    ph.run()
+
+    # Read forces and assemble the dynamical matrix
+    ph.read(acoustic=True)
+    ph.clean()
+
+    path = atoms.cell.bandpath('GXULGK', npoints=100)
+    bs = ph.get_band_structure(path)
+
+    dos = ph.get_dos(kpts=(20, 20, 20)).sample_grid(npts=100, width=1e-3)
+
+    forces = ph.get_force_constant()
+    print (forces)
+
+    # Plot the band structure and DOS:
+    import matplotlib.pyplot as plt
+    fig = plt.figure(1, figsize=(8, 4), dpi=300)
+    ax = fig.add_axes([.12, .07, .67, .85])
+
+    emax = 0.035
+
+    bs.plot(ax=ax, emin=-0.01, emax=emax)
+
+    dosax = fig.add_axes([.8, .07, .17, .85])
+    dosax.fill_between(dos.weights[0], dos.energy, y2=0, color='grey',
+                       edgecolor='k', lw=1)
+
+    dosax.set_ylim(-0.01, emax)
+    dosax.set_yticks([])
+    dosax.set_xticks([])
+    dosax.set_xlabel("DOS", fontsize=18)
+
+    fig.savefig('Al_phonon.png')
+    return
+
+interact_manual(calculate_phonons, x=FloatSlider(value=4.05, min=0, max=10, description='a (Angstrom)'));

TP-Vibrations/Si_DFT_gpaw_vibrations.py

@@ -0,0 +1,98 @@
+from gpaw import GPAW, PW
+from ase import Atoms
+from phonopy import Phonopy
+from phonopy.structure.atoms import PhonopyAtoms
+import numpy as np
+
+def get_gpaw(kpts_size=None):
+    if kpts_size is None:
+        calc = GPAW(mode=PW(300),
+                    kpts={'size': (2, 2, 2)})
+    else:
+        calc = GPAW(mode=PW(300),
+                    kpts={'size': kpts_size})
+    return calc
+
+def get_crystal():
+    a = 5.404
+    cell = PhonopyAtoms(symbols=(['Si'] * 8),
+                        cell=np.diag((a, a, a)),
+                        scaled_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)])
+    return cell
+
+def phonopy_pre_process(cell, supercell_matrix=None):
+
+    if supercell_matrix is None:
+        smat = [[2,0,0], [0,2,0], [0,0,2]],
+    else:
+        smat = supercell_matrix
+    phonon = Phonopy(cell,
+                     smat,
+                     primitive_matrix=[[0, 0.5, 0.5],
+                                       [0.5, 0, 0.5],
+                                       [0.5, 0.5, 0]])
+    phonon.generate_displacements(distance=0.03)
+    print("[Phonopy] Atomic displacements:")
+    disps = phonon.get_displacements()
+    for d in disps:
+        print("[Phonopy] %d %s" % (d[0], d[1:]))
+    return phonon
+
+def run_gpaw(calc, phonon):
+    supercells = phonon.get_supercells_with_displacements()
+    # Force calculations by calculator
+    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(("[Phonopy] Drift force:" + "%11.5f" * 3) % tuple(drift_force))
+        # Simple translational invariance
+        for force in forces:
+            force -= drift_force / forces.shape[0]
+        set_of_forces.append(forces)
+    return set_of_forces
+
+def phonopy_post_process(phonon, set_of_forces):
+    phonon.produce_force_constants(forces=set_of_forces)
+    print('')
+    print("[Phonopy] Phonon frequencies at Gamma:")
+    for i, freq in enumerate(phonon.get_frequencies((0, 0, 0))):
+        print("[Phonopy] %3d: %10.5f THz" %  (i + 1, freq)) # THz
+
+    # DOS
+    phonon.set_mesh([21, 21, 21])
+    phonon.set_total_DOS(tetrahedron_method=True)
+    print('')
+    print("[Phonopy] Phonon DOS:")
+    for omega, dos in np.array(phonon.get_total_DOS()).T:
+        print("%15.7f%15.7f" % (omega, dos))
+
+def main():
+    cell = get_crystal()
+
+    # 1x1x1 supercell of conventional unit cell
+    calc = get_gpaw(kpts_size=(4, 4, 4))
+    phonon = phonopy_pre_process(cell, supercell_matrix=np.eye(3, dtype='intc'))
+
+    # # 2x2x2 supercell of conventional unit cell
+    # calc = get_gpaw(kpts_size=(2, 2, 2))
+    # phonon = phonopy_pre_process(cell,
+    #                              supercell_matrix=(np.eye(3, dtype='intc') * 2))
+
+    set_of_forces = run_gpaw(calc, phonon)
+    phonopy_post_process(phonon, set_of_forces)
+
+if __name__ == "__main__":
+    main()