scnc
/
phonopy
mirror of https://github.com/phonopy/phonopy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'tdisp' into develop

This commit is contained in:
Atsushi Togo 2017-09-19 22:38:20 +09:00
parent 492805a932 5b96dfd9de
commit 1d2b532a89
11 changed files with 394 additions and 189 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
phonopy/api_phonopy.py
View File

@ -49,7 +49,7 @@ from phonopy.harmonic.dynamical_matrix import (DynamicalMatrix,
DynamicalMatrixNAC)
from phonopy.phonon.band_structure import BandStructure
from phonopy.phonon.thermal_properties import ThermalProperties
from phonopy.phonon.mesh import Mesh
from phonopy.phonon.mesh import Mesh, IterMesh
from phonopy.units import VaspToTHz
from phonopy.phonon.dos import TotalDos, PartialDos
from phonopy.phonon.thermal_displacement import (ThermalDisplacements,
@ -136,6 +136,9 @@ class Phonopy(object):
# set_mesh
self._mesh = None
# set_iter_mesh
self._iter_mesh = None
# set_tetrahedron_method
self._tetrahedron_method = None
@ -558,6 +561,56 @@ class Phonopy(object):
def write_yaml_mesh(self):
self._mesh.write_yaml()
# Sampling mesh:
# Solving dynamical matrices at q-points one-by-one as an iterator
def set_iter_mesh(self,
mesh,
shift=None,
is_time_reversal=True,
is_mesh_symmetry=True,
is_eigenvectors=False,
is_gamma_center=False):
"""Create an IterMesh instance
"""
if self._dynamical_matrix is None:
print("Warning: Dynamical matrix has not yet built.")
self._iter_mesh = None
return False
self._iter_mesh = IterMesh(
self._dynamical_matrix,
mesh,
shift=shift,
is_time_reversal=is_time_reversal,
is_mesh_symmetry=is_mesh_symmetry,
is_eigenvectors=is_eigenvectors,
is_gamma_center=is_gamma_center,
rotations=self._primitive_symmetry.get_pointgroup_operations(),
factor=self._factor)
return True
def get_iter_mesh(self):
"""Returns IterMesh instance
This instance object does not store all phonon data. With very
dense mesh and eigenvectors needed, IterMesh can save memory
space, but expected to be slow.
This object is used as an iterator. Phonon frequencies and
eigenvectos are obtained as below:
for i, (freqs, eigvecs) in enumerate(iter_mesh):
print(i + 1)
print(freqs)
print(eigvecs)
"""
return self._iter_mesh
# Plot band structure and DOS (PDOS) together
def plot_band_structure_and_dos(self, pdos_indices=None, labels=None):
import matplotlib.pyplot as plt
@ -799,37 +852,44 @@ class Phonopy(object):
are ignored.
direction:
Projection direction in reduced coordinates
Projection direction in reduced coordinates (
"""
self._thermal_displacements = None
if self._mesh is None:
print("Warning: \'set_mesh\' has to finish correctly "
"before \'set_thermal_displacements\'.")
return False
if self._mesh is not None:
eigvecs = self._mesh.get_eigenvectors()
mesh_nums = self._mesh.get_mesh_numbers()
if eigvecs is None:
print("Warning: Eigenvectors have to be calculated.")
return False
if np.prod(mesh_nums) != len(eigvecs):
print("Warning: Sampling mesh must not be symmetrized.")
return False
eigvecs = self._mesh.get_eigenvectors()
frequencies = self._mesh.get_frequencies()
mesh_nums = self._mesh.get_mesh_numbers()
iter_phonons = self._mesh
else:
if self._iter_mesh is not None:
iter_phonons = self._iter_mesh
else:
print("Warning: \'set_mesh\' has to finish correctly "
"before \'set_thermal_displacements\'.")
return False
if self._mesh.get_eigenvectors() is None:
print("Warning: Eigenvectors have to be calculated.")
return False
if direction is not None:
projection_direction = np.dot(direction, self._primitive.get_cell())
td = ThermalDisplacements(iter_phonons,
self._primitive.get_masses(),
projection_direction=projection_direction,
cutoff_frequency=cutoff_frequency)
else:
td = ThermalDisplacements(iter_phonons,
self._primitive.get_masses(),
cutoff_frequency=cutoff_frequency)
if np.prod(mesh_nums) != len(eigvecs):
print("Warning: Sampling mesh must not be symmetrized.")
return False
td = ThermalDisplacements(frequencies,
eigvecs,
self._primitive.get_masses(),
cutoff_frequency=cutoff_frequency)
if temperatures is None:
td.set_temperature_range(t_min, t_max, t_step)
else:
td.set_temperatures(temperatures)
if direction is not None:
td.project_eigenvectors(direction, self._primitive.get_cell())
td.run()
self._thermal_displacements = td
@ -875,28 +935,30 @@ class Phonopy(object):
"""
self._thermal_displacement_matrices = None
if self._mesh is None:
print("Warning: \'set_mesh\' has to finish correctly "
"before \'set_thermal_displacement_matrices\'.")
return False
if self._mesh is not None:
eigvecs = self._mesh.get_eigenvectors()
if eigvecs is None:
print("Warning: Eigenvectors have to be calculated.")
return False
if np.prod(self._mesh.get_mesh_numbers()) != len(eigvecs):
print("Warning: Sampling mesh must not be symmetrized.")
return False
eigvecs = self._mesh.get_eigenvectors()
frequencies = self._mesh.get_frequencies()
mesh_nums = self._mesh.get_mesh_numbers()
iter_phonons = self._mesh
else:
if self._iter_mesh is not None:
iter_phonons = self._iter_mesh
else:
print("Warning: \'set_mesh\' has to finish correctly "
"before \'set_thermal_displacement_matrices\'.")
return False
if self._mesh.get_eigenvectors() is None:
print("Warning: Eigenvectors have to be calculated.")
return False
tdm = ThermalDisplacementMatrices(
iter_phonons,
self._primitive.get_masses(),
cutoff_frequency=cutoff_frequency,
lattice=self._primitive.get_cell().T)
if np.prod(mesh_nums) != len(eigvecs):
print("Warning: Sampling mesh must not be symmetrized.")
return False
tdm = ThermalDisplacementMatrices(frequencies,
eigvecs,
self._primitive.get_masses(),
cutoff_frequency=cutoff_frequency,
lattice=self._primitive.get_cell().T)
if t_cif is None:
tdm.set_temperature_range(t_min, t_max, t_step)
else:
@ -1093,7 +1155,7 @@ class Phonopy(object):
return True
def get_modulated_supercells(self):
"""Returns cells with modulations as Atoms objects"""
"""Returns cells with modulations as Atoms instances"""
return self._modulation.get_modulated_supercells()
def get_modulations_and_supercell(self):
@ -1102,7 +1164,7 @@ class Phonopy(object):
(modulations, supercell)
modulations: Atomic modulations of supercell in Cartesian coordinates
supercell: Supercell as an Atoms object.
supercell: Supercell as an Atoms instance.
"""
return self._modulation.get_modulations_and_supercell()

4
phonopy/cui/settings.py
View File

@ -1563,7 +1563,7 @@ class PhonopyConfParser(ConfParser):
"TDISTANCE is incorrectly specified.")
if len(atom_pairs) > 0:
self.set_parameter('tdistance', atom_pairs)
# Projection direction used for thermal displacements and PDOS
if conf_key == 'projection_direction':
vals = [float(x) for x in confs['projection_direction'].split()]
@ -1829,7 +1829,7 @@ class PhonopyConfParser(ConfParser):
self._settings.set_projection_direction(
params['projection_direction'])
self._settings.set_is_mesh_symmetry(False)
# Thermal displacement matrices
if 'tdispmat' in params or 'tdispmat_cif' in params:
if 'tdispmat' in params and not params['tdispmat']:

129
phonopy/phonon/mesh.py
View File

@ -36,7 +36,7 @@ import numpy as np
from phonopy.units import VaspToTHz
from phonopy.structure.grid_points import GridPoints
class Mesh(object):
class MeshBase(object):
def __init__(self,
dynamical_matrix,
mesh,
@ -45,17 +45,13 @@ class Mesh(object):
is_mesh_symmetry=True,
is_eigenvectors=False,
is_gamma_center=False,
group_velocity=None,
rotations=None, # Point group operations in real space
factor=VaspToTHz,
use_lapack_solver=False):
factor=VaspToTHz):
self._mesh = np.array(mesh, dtype='intc')
self._is_eigenvectors = is_eigenvectors
self._factor = factor
self._cell = dynamical_matrix.get_primitive()
self._dynamical_matrix = dynamical_matrix
self._use_lapack_solver = use_lapack_solver
self._gp = GridPoints(self._mesh,
np.linalg.inv(self._cell.get_cell()),
@ -73,20 +69,12 @@ class Mesh(object):
self._eigenvalues = None
self._eigenvectors = None
self._group_velocity = group_velocity
self._group_velocities = None
def run(self):
self._set_phonon()
if self._group_velocity is not None:
self._set_group_velocities(self._group_velocity)
def get_dynamical_matrix(self):
return self._dynamical_matrix
def get_mesh_numbers(self):
return self._mesh
def get_qpoints(self):
return self._qpoints
@ -98,7 +86,7 @@ class Mesh(object):
def get_ir_grid_points(self):
return self._gp.get_ir_grid_points()
def get_grid_mapping_table(self):
return self._gp.get_grid_mapping_table()
@ -108,16 +96,69 @@ class Mesh(object):
def get_frequencies(self):
return self._frequencies
class Mesh(MeshBase):
def __init__(self,
dynamical_matrix,
mesh,
shift=None,
is_time_reversal=True,
is_mesh_symmetry=True,
is_eigenvectors=False,
is_gamma_center=False,
group_velocity=None,
rotations=None, # Point group operations in real space
factor=VaspToTHz,
use_lapack_solver=False):
MeshBase.__init__(self,
dynamical_matrix,
mesh,
shift=shift,
is_time_reversal=is_time_reversal,
is_mesh_symmetry=is_mesh_symmetry,
is_eigenvectors=is_eigenvectors,
is_gamma_center=is_gamma_center,
rotations=rotations,
factor=factor)
self._group_velocity = group_velocity
self._group_velocities = None
self._use_lapack_solver = use_lapack_solver
self._q_count = 0
def __iter__(self):
return self
def next(self):
return self.__next__()
def __next__(self):
if self._eigenvectors is None:
return StopIteration
if self._q_count == len(self._qpoints):
raise StopIteration
else:
i = self._q_count
self._q_count += 1
return self._frequencies[i], self._eigenvectors[i]
def run(self):
self._set_phonon()
if self._group_velocity is not None:
self._set_group_velocities(self._group_velocity)
def get_group_velocities(self):
return self._group_velocities
def get_eigenvectors(self):
"""
Eigenvectors is a numpy array of three dimension.
The first index runs through q-points.
In the second and third indices, eigenvectors obtained
using numpy.linalg.eigh are stored.
The third index corresponds to the eigenvalue's index.
The second index is for atoms [x1, y1, z1, x2, y2, z2, ...].
"""
@ -216,7 +257,55 @@ class Mesh(object):
dtype='double',
order='C') * self._factor
def _set_group_velocities(self, group_velocity):
group_velocity.set_q_points(self._qpoints)
self._group_velocities = group_velocity.get_group_velocity()
class IterMesh(MeshBase):
def __init__(self,
dynamical_matrix,
mesh,
shift=None,
is_time_reversal=True,
is_mesh_symmetry=True,
is_eigenvectors=False,
is_gamma_center=False,
rotations=None, # Point group operations in real space
factor=VaspToTHz):
MeshBase.__init__(self,
dynamical_matrix,
mesh,
shift=shift,
is_time_reversal=is_time_reversal,
is_mesh_symmetry=is_mesh_symmetry,
is_eigenvectors=is_eigenvectors,
is_gamma_center=is_gamma_center,
rotations=rotations,
factor=factor)
self._q_count = 0
def __iter__(self):
return self
def next(self):
return self.__next__()
def __next__(self):
if self._q_count == len(self._qpoints):
raise StopIteration
else:
q = self._qpoints[self._q_count]
self._dynamical_matrix.set_dynamical_matrix(q)
dm = self._dynamical_matrix.get_dynamical_matrix()
if self._is_eigenvectors:
eigvals, self._eigenvectors = np.linalg.eigh(dm)
self._eigenvalues = eigvals.real
else:
self._eigenvalues = np.linalg.eigvalsh(dm).real
self._frequencies = np.array(np.sqrt(abs(self._eigenvalues)) *
np.sign(self._eigenvalues),
dtype='double',
order='C') * self._factor
self._q_count += 1
return self._frequencies, self._eigenvectors

111
phonopy/phonon/thermal_displacement.py
View File

@ -40,8 +40,6 @@ from phonopy.interface.cif import write_cif_P1
class ThermalMotion(object):
def __init__(self,
frequencies, # have to be supplied in THz
eigenvectors,
masses,
cutoff_frequency=None):
@ -50,13 +48,8 @@ class ThermalMotion(object):
else:
self._cutoff_frequency = cutoff_frequency
self._distances = None
self._displacements = None
self._frequencies = frequencies
self._p_eigenvectors = None
self._eigenvectors = eigenvectors
self._masses = masses * AMU
self._masses3 = np.array([[m] * 3 for m in masses]).flatten() * AMU
self._masses3 = np.array([[m] * 3 for m in masses]).ravel() * AMU
self._temperatures = None
def get_Q2(self, freq, t): # freq in THz
@ -96,30 +89,6 @@ class ThermalMotion(object):
condition = np.logical_not(t_array < 0)
self._temperatures = np.extract(condition, t_array)
def project_eigenvectors(self, direction, lattice=None):
"""
direction
Without specifying lattice:
Projection direction in Cartesian coordinates
With lattice:
Projection direction in fractional coordinates
"""
if lattice is not None:
projector = np.dot(direction, lattice)
else:
projector = np.array(direction, dtype=float)
projector /= np.linalg.norm(projector)
self._p_eigenvectors = []
for vecs_q in self._eigenvectors:
p_vecs_q = []
for vecs in vecs_q.T:
p_vecs_q.append(np.dot(vecs.reshape(-1, 3), projector))
self._p_eigenvectors.append(np.transpose(p_vecs_q))
self._p_eigenvectors = np.array(self._p_eigenvectors)
def _get_population(self, freq, t): # freq in THz
if t < 1: # temperatue less than 1 K is approximated as 0 K.
return 0
@ -128,14 +97,29 @@ class ThermalMotion(object):
class ThermalDisplacements(ThermalMotion):
def __init__(self,
frequencies, # Have to be supplied in THz
eigenvectors,
iter_phonons,
masses,
projection_direction=None,
cutoff_frequency=None):
"""Calculate mean square displacements
iter_phonons: Mesh or IterMesh instance
masses: Atomic masses
projection_direction:
Eigenvector projection direction in Cartesian coordinates.
If None, eigenvector is not projected.
"""
self._iter_phonons = iter_phonons
if projection_direction is None:
self._projection_direction = None
else:
self._projection_direction = (projection_direction /
np.linalg.norm(projection_direction))
ThermalMotion.__init__(self,
frequencies,
eigenvectors,
masses,
cutoff_frequency=cutoff_frequency)
@ -145,25 +129,30 @@ class ThermalDisplacements(ThermalMotion):
return (self._temperatures, self._displacements)
def run(self):
freqs = self._frequencies
temps = self._temperatures
if self._p_eigenvectors is not None:
if self._projection_direction is not None:
masses = self._masses
eigvecs = self._p_eigenvectors
else:
masses = self._masses3
eigvecs = self._eigenvectors
temps = self._temperatures
disps = np.zeros((len(temps), len(masses)), dtype=float)
for fs, vecs2 in zip(freqs, abs(eigvecs) ** 2):
for f, v2 in zip(fs, vecs2.T):
for count, (fs, vecs) in enumerate(self._iter_phonons):
if self._projection_direction is not None:
p_vecs = []
for v in vecs.T:
p_vecs.append(np.dot(v.reshape(-1, 3),
self._projection_direction))
vecs2 = np.abs(p_vecs) ** 2
else:
vecs2 = (abs(vecs) ** 2).T
for f, v2 in zip(fs, vecs2):
if f > self._cutoff_frequency:
c = v2 / masses
for i, t in enumerate(temps):
disps[i] += self.get_Q2(f, t) * c
self._displacements = disps / len(freqs)
self._displacements = disps / (count + 1)
def write_yaml(self):
natom = len(self._masses)
f = open('thermal_displacements.yaml', 'w')
@ -192,22 +181,33 @@ class ThermalDisplacements(ThermalMotion):
if is_legend:
pyplot.legend(plots, labels, loc='upper left')
def _project_eigenvectors(self):
"""Eigenvectors are projected along Cartesian direction"""
self._p_eigenvectors = []
for vecs_q in self._eigenvectors:
p_vecs_q = []
for vecs in vecs_q.T:
p_vecs_q.append(np.dot(vecs.reshape(-1, 3),
self._projection_direction))
self._p_eigenvectors.append(np.transpose(p_vecs_q))
self._p_eigenvectors = np.array(self._p_eigenvectors)
class ThermalDisplacementMatrices(ThermalMotion):
def __init__(self,
frequencies, # Have to be supplied in THz
eigenvectors,
iter_phonons,
masses,
cutoff_frequency=None,
lattice=None): # column vectors in real space
ThermalMotion.__init__(self,
frequencies,
eigenvectors,
masses,
cutoff_frequency=cutoff_frequency)
self._disp_matrices = None
self._disp_matrices_cif = None
self._iter_phonons = iter_phonons
if lattice is not None:
A = lattice
@ -223,7 +223,7 @@ class ThermalDisplacementMatrices(ThermalMotion):
disps = np.zeros((len(self._temperatures), len(self._masses),
3, 3), dtype=complex)
for freqs, eigvecs in zip(self._frequencies, self._eigenvectors):
for count, (freqs, eigvecs) in enumerate(self._iter_phonons):
for f, vec in zip(freqs, eigvecs.T):
if f > self._cutoff_frequency:
c = []
@ -231,7 +231,7 @@ class ThermalDisplacementMatrices(ThermalMotion):
c.append(np.outer(v, v.conj()) / m)
for i, t in enumerate(self._temperatures):
disps[i] += self.get_Q2(f, t) * np.array(c)
self._disp_matrices = disps / len(self._frequencies)
self._disp_matrices = disps / (count + 1)
if self._ANinv is not None:
self._disp_matrices_cif = np.zeros(self._disp_matrices.shape,
@ -297,13 +297,16 @@ class ThermalDistances(ThermalMotion):
self._supercell = supercell
self._qpoints = qpoints
self._symprec = symprec
self._frequencies = frequencies
self._eigenvectors = eigenvectors
ThermalMotion.__init__(self,
frequencies,
eigenvectors,
primitive.get_masses(),
cutoff_frequency=cutoff_frequency)
self._p_eigenvectors = None
self._distances = None
def _get_cross(self, v, delta_r, q, atom1, atom2):
phase = np.exp(2j * np.pi * np.dot(delta_r, q))
cross_val = v[atom1]*phase*v[atom2].conjugate()

54
scripts/phonopy
View File

@ -731,30 +731,40 @@ elif run_mode == 'band' or run_mode == 'mesh' or run_mode == 'band_mesh':
t_symmetry,
q_symmetry,
is_gamma_center) = settings.get_mesh()
phonon.set_mesh(mesh,
mesh_shift,
is_time_reversal=t_symmetry,
is_mesh_symmetry=q_symmetry,
is_eigenvectors=settings.get_is_eigenvectors(),
is_gamma_center=settings.get_is_gamma_center(),
run_immediately=False)
weights = phonon.get_mesh()[1]
if log_level > 0:
if q_symmetry:
print("Number of irreducible q-points on sampling mesh: "
"%d/%d" % (weights.shape[0], np.prod(mesh)))
else:
print("Number of q-points on sampling mesh: %d" %
weights.shape[0])
print("Calculating phonons on sampling mesh...")
phonon.run_mesh()
if (settings.get_is_thermal_displacements() or
settings.get_is_thermal_displacement_matrices()):
phonon.set_iter_mesh(mesh,
mesh_shift,
is_time_reversal=t_symmetry,
is_mesh_symmetry=q_symmetry,
is_eigenvectors=settings.get_is_eigenvectors(),
is_gamma_center=settings.get_is_gamma_center())
else:
phonon.set_mesh(mesh,
mesh_shift,
is_time_reversal=t_symmetry,
is_mesh_symmetry=q_symmetry,
is_eigenvectors=settings.get_is_eigenvectors(),
is_gamma_center=settings.get_is_gamma_center(),
run_immediately=False)
weights = phonon.get_mesh()[1]
if log_level > 0:
if q_symmetry:
print("Number of irreducible q-points on sampling mesh: "
"%d/%d" % (weights.shape[0], np.prod(mesh)))
else:
print("Number of q-points on sampling mesh: %d" %
weights.shape[0])
print("Calculating phonons on sampling mesh...")
phonon.run_mesh()
if settings.get_write_mesh():
if settings.get_is_hdf5():
phonon.write_hdf5_mesh()
else:
phonon.write_yaml_mesh()
if settings.get_write_mesh():
if settings.get_is_hdf5():
phonon.write_hdf5_mesh()
else:
phonon.write_yaml_mesh()
# Thermal property
if settings.get_is_thermal_properties():

4
test/phonopy/BORN_NaCl Normal file
View File

@ -0,0 +1,4 @@
14.400
2.43533967 0 0 0 2.43533967 0 0 0 2.43533967
1.086875 0 0 0 1.086875 0 0 0 1.086875
-1.086875 0 0 0 -1.086875 0 0 0 -1.086875

0
test/phonopy/phonon/FORCE_SETS_moment → test/phonopy/FORCE_SETS_NaCl
View File

15
test/phonopy/phonon/POSCAR_moment
View File

@ -1,15 +0,0 @@
Na Cl
1.00000000000000
5.6903014761756712 0.0000000000000000 0.0000000000000000
0.0000000000000000 5.6903014761756712 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.6903014761756712
4 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000

10
test/phonopy/phonon/test_irreps.py
View File

@ -1,15 +1,16 @@
import unittest
import os
import sys
try:
from StringIO import StringIO
except ImportError:
from io import StringIO
import sys
import numpy as np
from phonopy import Phonopy
from phonopy.interface.vasp import read_vasp
from phonopy.file_IO import parse_FORCE_SETS
import os
data_dir=os.path.dirname(os.path.abspath(__file__))
data_dir = os.path.dirname(os.path.abspath(__file__))
chars_Pc = """ 1. 0. -1. 0.
1. 0. 1. 0.
@ -495,7 +496,8 @@ class TestIrreps(unittest.TestCase):
phonon = Phonopy(cell,
np.diag(dim),
primitive_matrix=pmat)
force_sets = parse_FORCE_SETS(filename=os.path.join(data_dir,"FORCE_SETS_%s" % spgtype))
filename = os.path.join(data_dir,"FORCE_SETS_%s" % spgtype)
force_sets = parse_FORCE_SETS(filename=filename)
phonon.set_displacement_dataset(force_sets)
phonon.produce_force_constants()
print(phonon.get_symmetry().get_pointgroup())

60
test/phonopy/phonon/test_mesh.py Normal file
View File

@ -0,0 +1,60 @@
import unittest
import os
import sys
try:
from StringIO import StringIO
except ImportError:
from io import StringIO
import numpy as np
from phonopy import Phonopy
from phonopy.interface.vasp import read_vasp
from phonopy.file_IO import parse_FORCE_SETS, parse_BORN
data_dir = os.path.dirname(os.path.abspath(__file__))
class TestIterMesh(unittest.TestCase):
def setUp(self):
pass
def tearDown(self):
pass
def testIterMesh(self):
phonon = self._get_phonon()
phonon.set_iter_mesh([3, 3, 3], is_eigenvectors=True)
imesh = phonon.get_iter_mesh()
freqs = []
eigvecs = []
for i, (f, e) in enumerate(imesh):
freqs.append(f)
eigvecs.append(e)
phonon.set_mesh([3, 3, 3], is_eigenvectors=True)
_, _, mesh_freqs, mesh_eigvecs = phonon.get_mesh()
np.testing.assert_allclose(mesh_freqs, freqs)
np.testing.assert_allclose(mesh_eigvecs, eigvecs)
def _get_phonon(self):
cell = read_vasp(os.path.join(data_dir, "../POSCAR_NaCl"))
phonon = Phonopy(cell,
np.diag([2, 2, 2]),
primitive_matrix=[[0, 0.5, 0.5],
[0.5, 0, 0.5],
[0.5, 0.5, 0]])
filename = os.path.join(data_dir,"../FORCE_SETS_NaCl")
force_sets = parse_FORCE_SETS(filename=filename)
phonon.set_displacement_dataset(force_sets)
phonon.produce_force_constants()
filename_born = os.path.join(data_dir, "../BORN_NaCl")
nac_params = parse_BORN(phonon.get_primitive(), filename=filename_born)
phonon.set_nac_params(nac_params)
return phonon
def _set_mesh(self):
pass
if __name__ == '__main__':
suite = unittest.TestLoader().loadTestsFromTestCase(TestIterMesh)
unittest.TextTestRunner(verbosity=2).run(suite)

48
test/phonopy/phonon/test_moment.py
View File

@ -8,22 +8,22 @@ import numpy as np
from phonopy import Phonopy
from phonopy.phonon.moment import PhononMoment
from phonopy.interface.vasp import read_vasp
from phonopy.file_IO import parse_FORCE_SETS
from phonopy.file_IO import parse_FORCE_SETS, parse_BORN
import os
data_dir=os.path.dirname(os.path.abspath(__file__))
result_full_range = """
1.000000 1.000000 1.000000
4.062877 4.237388 3.888351
17.935864 19.412878 16.458717
4.062876 4.237382 3.888355
17.935853 19.412818 16.458755
1.000000 1.000000 1.000000
3.515491 3.605997 3.436412
12.456606 13.099939 11.894498
3.515492 3.605997 3.436413
12.456609 13.099933 11.894505
"""
class TestMoment(unittest.TestCase):
def setUp(self):
self._cell = read_vasp(os.path.join(data_dir,"POSCAR_moment"))
self._cell = read_vasp(os.path.join(data_dir, "../POSCAR_NaCl"))
def tearDown(self):
pass
@ -52,6 +52,13 @@ class TestMoment(unittest.TestCase):
vals[i, 1:] = moms
self.assertTrue((np.abs(moms - data[i, 1:]) < 1e-5).all())
moment = PhononMoment(f, w)
moment.set_frequency_range(freq_min=3, freq_max=4)
for i in range(3):
moment.run(order=i)
vals[i + 3, 0] = moment.get_moment()
self.assertTrue(np.abs(moment.get_moment() - data[i + 3, 0]) < 1e-5)
moment = PhononMoment(f, w, eigenvectors=e)
moment.set_frequency_range(freq_min=3, freq_max=4)
for i in range(3):
@ -60,13 +67,6 @@ class TestMoment(unittest.TestCase):
vals[i + 3, 1:] = moms
self.assertTrue((np.abs(moms - data[i + 3, 1:]) < 1e-5).all())
moment = PhononMoment(f, w)
moment.set_frequency_range(freq_min=3, freq_max=4)
for i in range(3):
moment.run(order=i)
vals[i + 3, 0] = moment.get_moment()
self.assertTrue(np.abs(moment.get_moment() - data[i + 3, 0]) < 1e-5)
# self._show(vals)
def _show(self, vals):
@ -79,24 +79,14 @@ class TestMoment(unittest.TestCase):
primitive_matrix=[[0, 0.5, 0.5],
[0.5, 0, 0.5],
[0.5, 0.5, 0]])
force_sets = parse_FORCE_SETS(filename=os.path.join(data_dir,"FORCE_SETS_moment"))
filename = os.path.join(data_dir, "../FORCE_SETS_NaCl")
force_sets = parse_FORCE_SETS(filename=filename)
phonon.set_displacement_dataset(force_sets)
phonon.produce_force_constants()
supercell = phonon.get_supercell()
born_elems = {'Na': [[1.08703, 0, 0],
[0, 1.08703, 0],
[0, 0, 1.08703]],
'Cl': [[-1.08672, 0, 0],
[0, -1.08672, 0],
[0, 0, -1.08672]]}
born = [born_elems[s] for s in ['Na', 'Cl']]
epsilon = [[2.43533967, 0, 0],
[0, 2.43533967, 0],
[0, 0, 2.43533967]]
factors = 14.400
phonon.set_nac_params({'born': born,
'factor': factors,
'dielectric': epsilon})
filename_born = os.path.join(data_dir, "../BORN_NaCl")
nac_params = parse_BORN(phonon.get_primitive(), filename=filename_born)
phonon.set_nac_params(nac_params)
return phonon
if __name__ == '__main__':