Improve interface for getting parameters of NAC

This commit is contained in:
Atsushi Togo 2017-11-01 18:14:55 +09:00
parent 7e0fd74fc5
commit 1ad205d0f5
3 changed files with 95 additions and 99 deletions

View File

@ -42,9 +42,7 @@ import numpy as np
from phonopy.structure.atoms import PhonopyAtoms as Atoms
from phonopy.structure.atoms import symbol_map, atom_data
from phonopy.structure.cells import get_primitive, get_supercell
from phonopy.structure.symmetry import (Symmetry, get_site_symmetry,
get_pointgroup_operations)
from phonopy.harmonic.force_constants import similarity_transformation
from phonopy.structure.symmetry import Symmetry, symmetrize_borns_and_epsilon
from phonopy.file_IO import (write_FORCE_SETS, write_force_constants_to_hdf5,
write_FORCE_CONSTANTS)
@ -397,51 +395,6 @@ def get_born_OUTCAR(poscar_filename="POSCAR",
symmetrize_tensors=symmetrize_tensors,
symprec=symprec)
def symmetrize_borns_and_epsilon(borns,
epsilon,
ucell,
symprec=1e-5,
is_symmetry=True):
lattice = ucell.get_cell()
positions = ucell.get_scaled_positions()
u_sym = Symmetry(ucell, is_symmetry=is_symmetry, symprec=symprec)
rotations = u_sym.get_symmetry_operations()['rotations']
translations = u_sym.get_symmetry_operations()['translations']
ptg_ops = get_pointgroup_operations(rotations)
epsilon_ = symmetrize_2nd_rank_tensor(epsilon, ptg_ops, lattice)
for i, Z in enumerate(borns):
site_sym = get_site_symmetry(i,
lattice,
positions,
rotations,
translations,
symprec)
Z = symmetrize_2nd_rank_tensor(Z, site_sym, lattice)
borns_ = np.zeros_like(borns)
for i in range(len(borns)):
count = 0
for r, t in zip(rotations, translations):
count += 1
diff = np.dot(positions, r.T) + t - positions[i]
diff -= np.rint(diff)
dist = np.sqrt(np.sum(np.dot(diff, lattice) ** 2, axis=1))
j = np.nonzero(dist < symprec)[0][0]
r_cart = similarity_transformation(lattice.T, r)
borns_[i] += similarity_transformation(r_cart, borns[j])
borns_[i] /= count
return borns_, epsilon_
def symmetrize_2nd_rank_tensor(tensor, symmetry_operations, lattice):
sym_cart = [similarity_transformation(lattice.T, r)
for r in symmetry_operations]
sum_tensor = np.zeros_like(tensor)
for sym in sym_cart:
sum_tensor += similarity_transformation(sym, tensor)
return sum_tensor / len(symmetry_operations)
def _get_indep_borns(ucell,
borns,
epsilon,

View File

@ -36,6 +36,7 @@ import sys
import numpy as np
import phonopy.structure.spglib as spg
from phonopy.structure.atoms import PhonopyAtoms as Atoms
from phonopy.harmonic.force_constants import similarity_transformation
class Symmetry(object):
def __init__(self, cell, symprec=1e-5, is_symmetry=True):
@ -110,7 +111,7 @@ class Symmetry(object):
rotations = self._symmetry_operations['rotations']
translations = self._symmetry_operations['translations']
return get_site_symmetry(atom_number,
return self._get_site_symmetry(atom_number,
lattice,
positions,
rotations,
@ -129,6 +130,39 @@ class Symmetry(object):
"""
return self._reciprocal_operations
def _get_pointgroup_operations(self, rotations):
ptg_ops = []
for rot in rotations:
is_same = False
for tmp_rot in ptg_ops:
if (tmp_rot == rot).all():
is_same = True
break
if not is_same:
ptg_ops.append(rot)
return ptg_ops
def _get_site_symmetry(self,
atom_number,
lattice,
positions,
rotations,
translations,
symprec):
pos = positions[atom_number]
site_symmetries = []
for r, t in zip(rotations, translations):
rot_pos = np.dot(pos, r.T) + t
diff = pos - rot_pos
diff -= np.rint(diff)
diff = np.dot(diff, lattice)
if np.linalg.norm(diff) < symprec:
site_symmetries.append(r)
return np.array(site_symmetries, dtype='intc')
def _set_symmetry_dataset(self):
self._dataset = spg.get_symmetry_dataset(self._cell, self._symprec)
self._symmetry_operations = {
@ -180,7 +214,7 @@ class Symmetry(object):
def _set_pointgroup_operations(self):
rotations = self._symmetry_operations['rotations']
ptg_ops = get_pointgroup_operations(rotations)
ptg_ops = self._get_pointgroup_operations(rotations)
reciprocal_rotations = [rot.T for rot in ptg_ops]
exist_r_inv = False
for rot in ptg_ops:
@ -285,34 +319,42 @@ def get_lattice_vector_equivalence(point_symmetry):
return equivalence
def get_site_symmetry(atom_number,
lattice,
positions,
rotations,
translations,
symprec):
pos = positions[atom_number]
site_symmetries = []
def symmetrize_borns_and_epsilon(borns,
epsilon,
ucell,
symprec=1e-5,
is_symmetry=True):
lattice = ucell.get_cell()
positions = ucell.get_scaled_positions()
u_sym = Symmetry(ucell, is_symmetry=is_symmetry, symprec=symprec)
rotations = u_sym.get_symmetry_operations()['rotations']
translations = u_sym.get_symmetry_operations()['translations']
ptg_ops = u_sym.get_pointgroup_operations()
epsilon_ = _symmetrize_2nd_rank_tensor(epsilon, ptg_ops, lattice)
for i, Z in enumerate(borns):
site_sym = u_sym.get_site_symmetry(i)
Z = _symmetrize_2nd_rank_tensor(Z, site_sym, lattice)
borns_ = np.zeros_like(borns)
for i in range(len(borns)):
count = 0
for r, t in zip(rotations, translations):
rot_pos = np.dot(pos, r.T) + t
diff = pos - rot_pos
count += 1
diff = np.dot(positions, r.T) + t - positions[i]
diff -= np.rint(diff)
diff = np.dot(diff, lattice)
if np.linalg.norm(diff) < symprec:
site_symmetries.append(r)
dist = np.sqrt(np.sum(np.dot(diff, lattice) ** 2, axis=1))
j = np.nonzero(dist < symprec)[0][0]
r_cart = similarity_transformation(lattice.T, r)
borns_[i] += similarity_transformation(r_cart, borns[j])
borns_[i] /= count
return np.array(site_symmetries, dtype='intc')
return borns_, epsilon_
def get_pointgroup_operations(rotations):
ptg_ops = []
for rot in rotations:
is_same = False
for tmp_rot in ptg_ops:
if (tmp_rot == rot).all():
is_same = True
break
if not is_same:
ptg_ops.append(rot)
return ptg_ops
def _symmetrize_2nd_rank_tensor(tensor, symmetry_operations, lattice):
sym_cart = [similarity_transformation(lattice.T, r)
for r in symmetry_operations]
sum_tensor = np.zeros_like(tensor)
for sym in sym_cart:
sum_tensor += similarity_transformation(sym, tensor)
return sum_tensor / len(symmetry_operations)

View File

@ -51,8 +51,8 @@ parser = OptionParser()
parser.set_defaults(num_atoms=None,
primitive_axis =None,
supercell_matrix=None,
symmetrize_tensors=False,
read_vasprunxml=False,
symmetrize_tensors=True,
read_outcar=False,
symprec=1e-5)
parser.add_option("--dim", dest="supercell_matrix",
action="store", type="string",
@ -60,24 +60,25 @@ parser.add_option("--dim", dest="supercell_matrix",
parser.add_option("--pa", "--primitive_axis", dest="primitive_axis",
action="store", type="string",
help="Same as PRIMITIVE_AXIS tags")
parser.add_option("--st", "--symmetrize_tensors", dest="symmetrize_tensors",
action="store_true",
help="Symmetrize tensors")
parser.add_option("--nost", "--no_symmetrize_tensors",
dest="symmetrize_tensors",
action="store_false",
help="Prevent from symmetrizing tensors")
parser.add_option("--tolerance", dest="symprec", type="float",
help="Symmetry tolerance to search")
parser.add_option("--vasprunxml", dest="read_vasprunxml",
parser.add_option("--outcar", dest="read_outcar",
action="store_true",
help=("Read vasprun.xml instead of OUTCAR. "
"POSCAR is not necessary"))
help=("Read OUTCAR instead of vasprun.xml. "
"POSCAR is necessary in this case."))
(options, args) = parser.parse_args()
if args:
outcar_filename = args[0]
else:
if options.read_vasprunxml:
outcar_filename = "vasprun.xml"
else:
if options.read_outcar:
outcar_filename = "OUTCAR"
else:
outcar_filename = "vasprun.xml"
if len(args) > 1:
poscar_filename = args[1]
@ -111,17 +112,17 @@ with warnings.catch_warnings() as w:
warnings.simplefilter("always")
try:
if options.read_vasprunxml:
borns, epsilon, atom_indices = get_born_vasprunxml(
filename=outcar_filename,
if options.read_outcar:
borns, epsilon, atom_indices = get_born_OUTCAR(
poscar_filename=poscar_filename,
outcar_filename=outcar_filename,
primitive_matrix=primitive_axis,
supercell_matrix=supercell_matrix,
symmetrize_tensors=options.symmetrize_tensors,
symprec=options.symprec)
else:
borns, epsilon, atom_indices = get_born_OUTCAR(
poscar_filename=poscar_filename,
outcar_filename=outcar_filename,
borns, epsilon, atom_indices = get_born_vasprunxml(
filename=outcar_filename,
primitive_matrix=primitive_axis,
supercell_matrix=supercell_matrix,
symmetrize_tensors=options.symmetrize_tensors,