Python Phonopy.get_thermal_properties方法代碼示例

# 需要導入模塊: from phonopy import Phonopy [as 別名]
# 或者: from phonopy.Phonopy import get_thermal_properties [as 別名]
def get_phonopy_qha(energies, volumes, force_constants, structure, t_min, t_step, t_max, mesh, eos,
                      pressure=0):
    """
    Return phonopy QHA interface.

    Args:
        energies (list):
        volumes (list):
        force_constants (list):
        structure (Structure):
        t_min (float): min temperature
        t_step (float): temperature step
        t_max (float): max temperature
        mesh (list/tuple): reciprocal space density
        eos (str): equation of state used for fitting the energies and the volumes.
            options supported by phonopy: vinet, murnaghan, birch_murnaghan
        pressure (float): in GPa, optional.

    Returns:
        PhonopyQHA
    """
    from phonopy import Phonopy
    from phonopy.structure.atoms import Atoms as PhonopyAtoms
    from phonopy import PhonopyQHA
    from phonopy.units import EVAngstromToGPa

    phon_atoms = PhonopyAtoms(symbols=[str(s.specie) for s in structure],
                              scaled_positions=structure.frac_coords,
                              cell=structure.lattice.matrix)
    scell = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    phonon = Phonopy(phon_atoms, scell)
    # compute the required phonon thermal properties
    temperatures = []
    free_energy = []
    entropy = []
    cv = []
    for f in force_constants:
        phonon.set_force_constants(-np.array(f))
        phonon.set_mesh(list(mesh))
        phonon.set_thermal_properties(t_step=t_step, t_min=t_min, t_max=t_max)
        t, g, e, c = phonon.get_thermal_properties()
        temperatures.append(t)
        free_energy.append(g)
        entropy.append(e)
        cv.append(c)

    # add pressure contribution
    energies = np.array(energies) + np.array(volumes) * pressure / EVAngstromToGPa
    # quasi-harmonic approx
    return PhonopyQHA(volumes, energies, eos=eos, temperatures=temperatures[0],
                      free_energy=np.array(free_energy).T, cv=np.array(cv).T,
                      entropy=np.array(entropy).T, t_max=np.max(temperatures[0]))

# 需要導入模塊: from phonopy import Phonopy [as 別名]
# 或者: from phonopy.Phonopy import get_thermal_properties [as 別名]
for q, d, freq in zip(q_points, distances, frequencies):
    print q, d, freq
phonon.plot_band_structure().show()

# Mesh sampling 20x20x20
phonon.set_mesh([20, 20, 20])
phonon.set_thermal_properties(t_step=10,
                              t_max=1000,
                              t_min=0)

# DOS
phonon.set_total_DOS(sigma=0.1)
for omega, dos in np.array(phonon.get_total_DOS()).T:
    print "%15.7f%15.7f" % (omega, dos)
phonon.plot_total_DOS().show()

# Thermal properties
for t, free_energy, entropy, cv in np.array(phonon.get_thermal_properties()).T:
    print ("%12.3f " + "%15.7f" * 3) % ( t, free_energy, entropy, cv )
phonon.plot_thermal_properties().show()

# PDOS
phonon.set_mesh([10, 10, 10],
                is_mesh_symmetry=False,
                is_eigenvectors=True)
phonon.set_partial_DOS(tetrahedron_method=True)
omegas, pdos = phonon.get_partial_DOS()
pdos_indices = [[0], [1]]
phonon.plot_partial_DOS(pdos_indices=pdos_indices,
                        legend=pdos_indices).show()

# 需要導入模塊: from phonopy import Phonopy [as 別名]
# 或者: from phonopy.Phonopy import get_thermal_properties [as 別名]

#.........這裏部分代碼省略.........
         phonon.set_thermal_properties(t_step=10,
                                       t_max=2500,
                                       t_min=0)
     
     
     elif species == 'Re': 
         fc = vasp.get_force_constants_vasprun_xml(vasprun,9,0)
         s = 5.
         a = superc.get_cell()[0][0]*2.
         print a
         bulk = PhonopyAtoms(symbols=['Re'] * 1,
                             scaled_positions= primitive.get_scaled_positions())
         bulk.set_cell(np.diag((a, a, a)))
         phonon = Phonopy(bulk,
                          [[s,0.,0.],[0.,s,0.],[0.,0.,s]],
                          primitive_matrix=[[-0.5, 0.5, 0.5],[0.5, -0.5, 0.5],[0.5, 0.5, -0.5]],
                          distance=0.01, factor=15.633302)
         print fc
         phonon.set_force_constants(fc[0])
         phonon.set_dynamical_matrix()
         #print phonon.get_dynamical_matrix_at_q([0,0,0])
         mesh = [100, 100, 100]
         phonon.set_mesh(mesh)
         qpoints, weights, frequencies, eigvecs = phonon.get_mesh()
         print frequencies
         phonon.set_total_DOS()
         
         phonon.set_thermal_properties(t_step=10,
                                       t_max=2500,
                                       t_min=0)
     
     
     f = open('F_TV','w')
     for t, free_energy, entropy, cv in np.array(phonon.get_thermal_properties()).T:
         #print t, cv
         #print ("%12.3f " + "%15.7f" * 3) % ( t, free_energy, entropy, cv )
         f.write(("%12.3f " + "%15.7f" + "\n") % ( t, free_energy))
     f.close()
     
     fc = open('thermal_properties','w')
     for t, free_energy, entropy, cv in np.array(phonon.get_thermal_properties()).T:
         fc.write(("%12.3f " + "%15.7f" *3 + "\n") % ( t, free_energy, entropy, cv ))
     fc.close()
     
     #phonon.plot_thermal_properties().show()
     
     #phonon.plot_total_DOS().show()
     phonon.write_total_DOS()
     #phonon.write_partial_DOS()
     phonon.write_yaml_thermal_properties()
     
     bands = []
     
     #### PRIMITIVE
     
     q_start  = np.array([0.0, 0.0, 0.0])
     #q_start  = np.array([0.5, 0.5, 0.0])
     q_end    = np.array([-0.5, 0.5, 0.5])
     #q_end    = np.array([0., 0., 0.])
     band = []
     for i in range(101):
         band.append(q_start + (q_end - q_start) / 100 * i)
     bands.append(band)
     
     band = []
     

# 需要導入模塊: from phonopy import Phonopy [as 別名]
# 或者: from phonopy.Phonopy import get_thermal_properties [as 別名]
print "%12s %15s%15s%15s" % ('T [K]',
                              'F [kJ/mol]',
                              'S [J/K/mol]',
                              'C_v [J/K/mol]')

# get_thermal_properties returns numpy array of
#
# [[ temperature, free energy, entropy, heat capacity ],
#  [ temperature, free energy, entropy, heat capacity ],...,]
#
# Frequency has to be given in THz internally. Therefore unit
# conversion factor may be specified when calling Phonon class. The
# unit of frequency in a calculator is square root of the unit of
# dynamical matrix, i.e.,
#
# /      [energy]       \^(1/2)
# | ------------------- |
# \ [mass] [distance]^2 /
#
# THz is the value above divided by 2pi*1e12 (2pi comes from the
# factor between angular frequency and ordinary frequency). See
# units.py in the phonopy directory.
phonon.set_mesh( mesh, shift )
phonon.set_thermal_properties( t_step=10,
                               t_max=1000,
                               t_min=0 )
for t, free_energy, entropy, cv in phonon.get_thermal_properties():
    print ("%12.3f " + "%15.7f" * 3) % ( t, free_energy, entropy, cv )

phonon.plot_thermal_properties().show()