diff --git a/phonopy/api_phonopy.py b/phonopy/api_phonopy.py index 47f316ca..c9560b36 100644 --- a/phonopy/api_phonopy.py +++ b/phonopy/api_phonopy.py @@ -37,7 +37,7 @@ import warnings import numpy as np from phonopy.version import __version__ from phonopy.interface import PhonopyYaml -from phonopy.structure.atoms import PhonopyAtoms as Atoms +from phonopy.structure.atoms import PhonopyAtoms from phonopy.structure.symmetry import Symmetry, symmetrize_borns_and_epsilon from phonopy.structure.cells import get_supercell, get_primitive from phonopy.harmonic.displacement import (get_least_displacements, @@ -104,7 +104,7 @@ class Phonopy(object): self._log_level = log_level # Create supercell and primitive cell - self._unitcell = Atoms(atoms=unitcell) + self._unitcell = PhonopyAtoms(atoms=unitcell) self._supercell_matrix = supercell_matrix self._primitive_matrix = primitive_matrix self._supercell = None @@ -250,6 +250,27 @@ class Phonopy(object): @property def displacements(self): + """Return displacements + + Returns + ------- + There are two types of displacement dataset. See the docstring + of set_displacement_dataset about types 1 and 2 for displacement + dataset format. + + Type-1, List of list + The internal list has 4 elements such as [32, 0.01, 0.0, 0.0]]. + The first element is the supercell atom index starting with 0. + The remaining three elements give the displacement in Cartesian + coordinates. + Type-2, array_like + Displacements of all atoms of all supercells in Cartesian + coordinates. + shape=(supercells, natom, 3) + dtype='double' + + """ + disps = [] if 'first_atoms' in self._displacement_dataset: for disp in self._displacement_dataset['first_atoms']: @@ -270,9 +291,9 @@ class Phonopy(object): Parameters ---------- displacemens : array_like - Snapshots of atomic displacements of all atoms in supercell. + Atomic displacements of all atoms of all supercells. Only all displacements in each supercell case is supported. - shape=(snapshots, natom, 3) + shape=(supercells, natom, 3) dtype='double' order='C' @@ -303,10 +324,10 @@ class Phonopy(object): forces = [] for disp in self._displacement_dataset['first_atoms']: if 'forces' in disp: - forces.append = disp['forces'] + forces.append(disp['forces']) return forces else: - return [] + return None @property def dynamical_matrix(self): @@ -429,22 +450,36 @@ class Phonopy(object): ---------- displacement_dataset : dict There are two dict structures. - 1. One atomic displacement in each supercell: + Type 1. One atomic displacement in each supercell: {'natom': number of atoms in supercell, 'first_atoms': [ {'number': atom index of displaced atom, 'displacement': displacement in Cartesian coordinates, 'forces': forces on atoms in supercell}, {...}, ...]} - 2. All atomic displacements in each supercell: + Type 2. All atomic displacements in each supercell: {'natom': number of atoms in supercell, 'displacements': ndarray, dtype='double', order='C', - shape=(natom, snapshots, 3) + shape=(supercells, natom, 3) 'forces': ndarray, dtype='double',, order='C', - shape=(natom, snapshots, 3)} + shape=(supercells, natom, 3)} + In type 2, displacements and forces can be given by numpy array + with different shape but that can be reshaped to + (supercells, natom, 3). """ - self._displacement_dataset = displacement_dataset + dds = displacement_dataset + if 'displacements' in dds: + natom = self._supercell.get_number_of_atoms() + if type(dds['displacements']) is np.ndarray: + if dds['displacements'].ndim in (1, 2): + d = dds['displacements'].reshape((-1, natom, 3)) + dds['displacements'] = d + if type(dds['forces']) is np.ndarray: + if dds['forces'].ndim in (1, 2): + f = dds['forces'].reshape((-1, natom, 3)) + dds['forces'] = f + self._displacement_dataset = dds self._supercells_with_displacements = None @forces.setter @@ -1527,7 +1562,7 @@ class Phonopy(object): return True def get_modulated_supercells(self): - """Returns cells with modulations as Atoms instances""" + """Returns cells with modulations as PhonopyAtoms instances""" return self._modulation.get_modulated_supercells() def get_modulations_and_supercell(self): @@ -1536,7 +1571,7 @@ class Phonopy(object): (modulations, supercell) modulations: Atomic modulations of supercell in Cartesian coordinates - supercell: Supercell as an Atoms instance. + supercell: Supercell as an PhonopyAtoms instance. """ return self._modulation.get_modulations_and_supercell() @@ -1749,7 +1784,7 @@ class Phonopy(object): for disp in self._displacement_dataset['first_atoms']: positions = self._supercell.get_positions() positions[disp['number']] += disp['displacement'] - supercells.append(Atoms( + supercells.append(PhonopyAtoms( numbers=self._supercell.get_atomic_numbers(), masses=self._supercell.get_masses(), magmoms=self._supercell.get_magnetic_moments(), diff --git a/phonopy/file_IO.py b/phonopy/file_IO.py index 3789e203..22efbda4 100644 --- a/phonopy/file_IO.py +++ b/phonopy/file_IO.py @@ -50,6 +50,22 @@ def write_FORCE_SETS(dataset, filename='FORCE_SETS'): def get_FORCE_SETS_lines(dataset, forces=None): + """Generate FORCE_SETS string + + See the format of dataset in the docstring of + Phonopy.set_displacement_dataset. Optionally for the type-1 (traditional) + format, forces can be given. In this case, sets of forces are + unnecessary to be stored in the dataset. + + """ + + if 'first_atoms' in dataset: + return _get_FORCE_SETS_lines_type1(dataset, forces=forces) + elif 'forces' in dataset: + return _get_FORCE_SETS_lines_type2(dataset) + + +def _get_FORCE_SETS_lines_type1(dataset, forces=None): num_atom = dataset['natom'] displacements = dataset['first_atoms'] if forces is None: @@ -69,6 +85,17 @@ def get_FORCE_SETS_lines(dataset, forces=None): return lines + +def _get_FORCE_SETS_lines_type2(dataset): + lines = [] + for displacements, forces in zip(dataset['displacements'], + dataset['forces']): + for d, f in zip(displacements, forces): + lines.append(("%15.8f" * 6) % (tuple(d) + tuple(f))) + + return lines + + def parse_FORCE_SETS(natom=None, is_translational_invariance=False, filename="FORCE_SETS"): @@ -103,14 +130,19 @@ def _get_set_of_forces(f, natom=None, is_translational_invariance=False): def _get_set_of_forces_type2(f, natom): data = np.loadtxt(f, dtype='double') - if data.shape[1] != 6 or data.shape[0] % natom != 0: + if data.shape[1] != 6 or (natom and data.shape[0] % natom != 0): msg = "Data shape of forces and displacements is incorrect." raise RuntimeError(msg) - data = data.reshape(-1, natom, 6) - dataset = {'natom': natom} - dataset['displacements'] = np.array(data[:, :, :3], - dtype='double', order='C') - dataset['forces'] = np.array(data[:, :, 3:], dtype='double', order='C') + if natom: + data = data.reshape(-1, natom, 6) + displacements = data[:, :, :3] + forces = data[:, :, 3:] + else: + displacements = data[:, :3] + forces = data[:, 3:] + dataset = {'displacements': + np.array(displacements, dtype='double', order='C'), + 'forces': np.array(forces, dtype='double', order='C')} return dataset diff --git a/phonopy/harmonic/force_constants.py b/phonopy/harmonic/force_constants.py index 009f8f6b..6423dcdd 100644 --- a/phonopy/harmonic/force_constants.py +++ b/phonopy/harmonic/force_constants.py @@ -284,7 +284,7 @@ def get_rotated_displacement(displacements, site_sym_cart): rot_disps = [] for u in displacements: rot_disps.append([np.dot(sym, u) for sym in site_sym_cart]) - return np.reshape(rot_disps, (-1, 3)) + return np.array(np.reshape(rot_disps, (-1, 3)), dtype='double', order='C') def get_rotated_forces(forces_syms, site_sym_cart): diff --git a/phonopy/interface/__init__.py b/phonopy/interface/__init__.py index 23678303..26fb58a6 100644 --- a/phonopy/interface/__init__.py +++ b/phonopy/interface/__init__.py @@ -48,13 +48,14 @@ def get_interface_mode(args): """ calculator_list = ['wien2k', 'abinit', 'qe', 'elk', 'siesta', 'cp2k', - 'crystal', 'vasp','dftbp'] + 'crystal', 'vasp', 'dftbp'] for calculator in calculator_list: mode = "%s_mode" % calculator if mode in args and args.__dict__[mode]: return calculator return None + def write_supercells_with_displacements(interface_mode, supercell, cells_with_disps, @@ -188,6 +189,7 @@ def read_crystal_structure(filename=None, unitcell = read_dftbp(cell_filename) return unitcell, (cell_filename,) + def get_default_cell_filename(interface_mode): if interface_mode == 'phonopy_yaml': return "phonopy_disp.yaml" @@ -210,6 +212,7 @@ def get_default_cell_filename(interface_mode): else: return None + def get_default_supercell_filename(interface_mode): if interface_mode == 'phonopy_yaml': return "phonopy_disp.yaml" @@ -226,7 +229,7 @@ def get_default_supercell_filename(interface_mode): elif interface_mode == 'crystal': return None # supercell.ext can not be parsed by crystal interface. elif interface_mode == 'dftbp': - return "geo.gen" + return "geo.genS" else: return None @@ -319,7 +322,7 @@ def get_default_physical_units(interface_mode): units['nac_factor'] = Hartree * Bohr units['distance_to_A'] = Bohr units['force_constants_unit'] = 'hartree/au^2' - units['length_unit'] = 'Angstrom' + units['length_unit'] = 'au' return units diff --git a/phonopy/interface/phonopy_yaml.py b/phonopy/interface/phonopy_yaml.py index 4d791c8c..b054941b 100644 --- a/phonopy/interface/phonopy_yaml.py +++ b/phonopy/interface/phonopy_yaml.py @@ -248,15 +248,37 @@ class PhonopyYaml(object): return lines def _force_sets_yaml_lines(self): + if 'first_atoms' in self._dataset: + return self._force_sets_yaml_lines_type1() + elif 'displacements' in self._dataset: + return self._force_sets_yaml_lines_type2() + else: + return [] + + def _force_sets_yaml_lines_type1(self): lines = ["displacements:", ] for i, d in enumerate(self._dataset['first_atoms']): lines.append("- atom: %4d" % (d['number'] + 1)) lines.append(" displacement:") lines.append(" [ %20.16f,%20.16f,%20.16f ]" % tuple(d['displacement'])) - lines.append(" forces:") - for f in d['forces']: - lines.append(" - [ %20.16f,%20.16f,%20.16f ]" % tuple(f)) + if 'forces' in d: + lines.append(" forces:") + for f in d['forces']: + lines.append(" - [ %20.16f,%20.16f,%20.16f ]" % tuple(f)) + lines.append("") + return lines + + def _force_sets_yaml_lines_type2(self): + lines = ["displacements:", ] + for i, (dset, fset) in enumerate(zip(self._dataset['displacements'], + self._dataset['forces'])): + lines.append("- # %4d" % (i + 1)) + for j, (d, f) in enumerate(zip(dset, fset)): + lines.append(" - displacement: # %d" % (j + 1)) + lines.append(" [ %20.16f,%20.16f,%20.16f ]" % tuple(d)) + lines.append(" force:") + lines.append(" [ %20.16f,%20.16f,%20.16f ]" % tuple(f)) lines.append("") return lines @@ -294,8 +316,11 @@ class PhonopyYaml(object): ('points' in self.yaml or 'atoms' in self.yaml)): self.unitcell = self._parse_cell(self.yaml) if 'displacements' in self.yaml: - if 'forces' in self.yaml['displacements'][0]: - self.dataset = self._parse_force_sets() + disp = self.yaml['displacements'][0] + if type(disp) is dict and 'forces' in disp: + self.dataset = self._parse_force_sets_type1() + elif type(disp) is list and 'forces' in disp[0]: + self.dataset = self._parse_force_sets_type2() if 'supercell_matrix' in self.yaml: self.supercell_matrix = np.array(self.yaml['supercell_matrix'], dtype='intc', order='C') @@ -368,7 +393,7 @@ class PhonopyYaml(object): else: return None - def _parse_force_sets(self): + def _parse_force_sets_type1(self): natom = len(self.yaml['displacements'][0]['forces']) dataset = {'natom': natom} first_atoms = [] @@ -379,3 +404,14 @@ class PhonopyYaml(object): 'forces': np.array(d['forces'], dtype='double', order='C')}) dataset['first_atoms'] = first_atoms return dataset + + def _parse_force_sets_type2(self): + nsets = len(self.yaml['displacements']) + natom = len(self.yaml['displacements'][0]) + forces = np.zeros((nsets, natom, 3), dtype='double', order='C') + displacements = np.zeros((nsets, natom, 3), dtype='double', order='C') + for i, dfset in enumerate(self.yaml['displacements']): + for j, df in enumerate(dfset): + forces[i, j] = df['force'] + displacements[i, j] = df['displacement'] + return {'forces': forces, 'displacements': displacements} diff --git a/phonopy/structure/atoms.py b/phonopy/structure/atoms.py index c10c7d82..14005dd1 100644 --- a/phonopy/structure/atoms.py +++ b/phonopy/structure/atoms.py @@ -36,11 +36,13 @@ import warnings import numpy as np + def Atoms(*args, **kwargs): warnings.warn("phonopy.atoms.Atoms is deprecated. Please use " "PhonopyAtoms instead of Atoms.", DeprecationWarning) return PhonopyAtoms(*args, **kwargs) + class _Atoms(object): """A class compatible with the ASE Atoms class. @@ -266,8 +268,18 @@ class PhonopyAtoms(_Atoms): symbols=self.symbols, pbc=True) + def totuple(self): + if self.magmoms: + return (self.cell, self.scaled_positions, self.numbers, + self.magmoms) + else: + return (self.cell, self.scaled_positions, self.numbers) + def to_tuple(self): - return (self.cell, self.scaled_positions, self.numbers) + warnings.warn( + "PhonopyAtoms.to_tuple is deprecated. Please use " + "PhonopyAtoms.totuple instead.", DeprecationWarning) + return self.totuple() def get_yaml_lines(self): lines = ["lattice:"] diff --git a/phonopy/structure/symmetry.py b/phonopy/structure/symmetry.py index 79d6688c..8241cad1 100644 --- a/phonopy/structure/symmetry.py +++ b/phonopy/structure/symmetry.py @@ -181,7 +181,8 @@ class Symmetry(object): return np.array(site_symmetries, dtype='intc') def _set_symmetry_dataset(self): - self._dataset = spg.get_symmetry_dataset(self._cell, self._symprec) + self._dataset = spg.get_symmetry_dataset(self._cell.totuple(), + self._symprec) self._symmetry_operations = { 'rotations': self._dataset['rotations'], 'translations': self._dataset['translations']} @@ -192,11 +193,7 @@ class Symmetry(object): self._map_atoms = self._dataset['equivalent_atoms'] def _set_symmetry_operations_with_magmoms(self): - cell = (self._cell.get_cell(), - self._cell.get_scaled_positions(), - self._cell.get_atomic_numbers(), - self._cell.get_magnetic_moments()) - self._symmetry_operations = spg.get_symmetry(cell, + self._symmetry_operations = spg.get_symmetry(self._cell.totuple(), symprec=self._symprec) self._map_atoms = self._symmetry_operations['equivalent_atoms'] self._set_map_atoms() @@ -302,7 +299,7 @@ def find_primitive(cell, symprec=1e-5): cell is found, an object of Atoms class of the primitive cell is returned. When not, None is returned. """ - lattice, positions, numbers = spg.find_primitive(cell, symprec) + lattice, positions, numbers = spg.find_primitive(cell.totuple(), symprec) if lattice is None: return None else: diff --git a/scripts/phonopy b/scripts/phonopy index d6bd0684..3f5d3763 100755 --- a/scripts/phonopy +++ b/scripts/phonopy @@ -555,7 +555,11 @@ else: # Read FORCE_SETS, FORCE_CONSTANTS, or force_constants.hdf5 else: file_exists("FORCE_SETS", log_level) force_sets = parse_FORCE_SETS(natom=num_satom) - if force_sets['natom'] != num_satom: + if 'natom' in force_sets: + natom = force_sets['natom'] + else: + natom = force_sets['forces'].shape[1] + if natom != num_satom: error_text = "Number of atoms in supercell is not consistent with " error_text += "the data in FORCE_SETS.\n" error_text += ("Please carefully check DIM, FORCE_SETS,"