Python symmetry.Symmetry类代码示例

    def __init__(self,
                 fc4,
                 supercell,
                 primitive,
                 mesh,
                 temperatures=None,
                 band_indices=None,
                 frequency_factor_to_THz=VaspToTHz,
                 is_nosym=False,
                 symprec=1e-3,
                 cutoff_frequency=1e-4,
                 log_level=False,
                 lapack_zheev_uplo='L'):
        self._fc4 = fc4
        self._supercell = supercell
        self._primitive = primitive
        self._masses = np.double(self._primitive.get_masses())
        self._mesh = np.intc(mesh)
        if temperatures is None:
            self._temperatures = np.double([0])
        else:
            self._temperatures = np.double(temperatures)
        num_band = primitive.get_number_of_atoms() * 3
        if band_indices is None:
            self._band_indices = np.arange(num_band, dtype='intc')
        else:
            self._band_indices = np.intc(band_indices)
        self._frequency_factor_to_THz = frequency_factor_to_THz
        self._is_nosym = is_nosym
        self._symprec = symprec
        self._cutoff_frequency = cutoff_frequency
        self._log_level = log_level
        self._lapack_zheev_uplo = lapack_zheev_uplo

        symmetry = Symmetry(primitive, symprec=symprec)
        self._point_group_operations = symmetry.get_pointgroup_operations()

        self._grid_address = None
        self._bz_map = None
        self._set_grid_address()
        
        self._grid_point = None
        self._quartets_at_q = None
        self._weights_at_q = None

        self._phonon_done = None
        self._frequencies = None
        self._eigenvectors = None
        self._dm = None
        self._nac_q_direction = None

        self._frequency_shifts = None
        
        # Unit to THz of Delta
        self._unit_conversion = (EV / Angstrom ** 4 / AMU ** 2
                                 / (2 * np.pi * THz) ** 2
                                 * Hbar * EV / (2 * np.pi * THz) / 8
                                 / np.prod(self._mesh))

        self._allocate_phonon()

def read_crystal(filename):
    f_crystal = open(filename)
    crystal_in = CrystalIn(f_crystal.readlines())
    f_crystal.close()
    tags = crystal_in.get_tags()
     
     
    cell = Atoms(cell=tags['lattice_vectors'], 
                 symbols=tags['atomic_species'],
                 scaled_positions=tags['coordinates'])
    
    magmoms = tags['magnetic_moments']
    if magmoms is not None:
        # Print out symmetry information for magnetic cases
        # Original code from structure/symmetry.py
        symmetry = Symmetry(cell, symprec=1e-5)
        print("CRYSTAL-interface: Magnetic structure, number of operations without spin: %d" %
              len(symmetry.get_symmetry_operations()['rotations']))
        print("CRYSTAL-interface: Spacegroup without spin: %s" % symmetry.get_international_table())

        cell.set_magnetic_moments(magmoms)
        symmetry = Symmetry(cell, symprec=1e-5)
        print("CRYSTAL-interface: Magnetic structure, number of operations with spin: %d" %
              len(symmetry.get_symmetry_operations()['rotations']))
        print("")
    
    return cell, tags['conv_numbers']

def get_born_OUTCAR(poscar_filename="POSCAR",
                    outcar_filename="OUTCAR",
                    primitive_axis=np.eye(3),
                    is_symmetry=True,
                    symmetrize_tensors=False):
    cell = read_vasp(poscar_filename)
    primitive = Primitive(cell, primitive_axis)
    p2p = primitive.get_primitive_to_primitive_map()
    symmetry = Symmetry(primitive, is_symmetry=is_symmetry)
    independent_atoms = symmetry.get_independent_atoms()
    prim_lat = primitive.get_cell().T
    outcar = open(outcar_filename)
    
    borns = []
    while True:
        line = outcar.readline()
        if not line:
            break
    
        if "NIONS" in line:
            num_atom = int(line.split()[11])
    
        if "MACROSCOPIC STATIC DIELECTRIC TENSOR" in line:
            epsilon = []
            outcar.readline()
            epsilon.append([float(x) for x in outcar.readline().split()])
            epsilon.append([float(x) for x in outcar.readline().split()])
            epsilon.append([float(x) for x in outcar.readline().split()])
    
        if "BORN" in line:
            outcar.readline()
            line = outcar.readline()
            if "ion" in line:
                for i in range(num_atom):
                    born = []
                    born.append([float(x) for x in outcar.readline().split()][1:])
                    born.append([float(x) for x in outcar.readline().split()][1:])
                    born.append([float(x) for x in outcar.readline().split()][1:])
                    outcar.readline()
                    borns.append(born)


    reduced_borns = []
    for p_i, u_i in enumerate(p2p):
        if p_i in independent_atoms:
            if symmetrize_tensors:
                site_sym = [similarity_transformation(prim_lat, rot)
                            for rot in symmetry.get_site_symmetry(p_i)]
                reduced_borns.append(symmetrize_tensor(borns[u_i], site_sym))
            else:
                reduced_borns.append(borns[u_i])
                
    if symmetrize_tensors:
        point_sym = [similarity_transformation(prim_lat, rot)
                     for rot in symmetry.get_pointgroup_operations()]
        epsilon = symmetrize_tensor(epsilon, point_sym)
    else:
        epsilon = np.array(epsilon)

    return np.array(reduced_borns), epsilon

def write_supercells_with_displacements( supercell,
                                         cells_with_displacements,
                                         npts, r0s, rmts,
                                         supercell_matrix,
                                         filename="wien2k-" ):
    v = supercell_matrix
    det = v[0,0] * v[1,1] * v[2,2] \
        + v[0,1] * v[1,2] * v[2,0] \
        + v[0,2] * v[1,0] * v[2,1] \
        - v[0,0] * v[1,2] * v[2,1] \
        - v[0,1] * v[1,0] * v[2,2] \
        - v[0,2] * v[1,1] * v[2,0]
        
    npts_super = []
    r0s_super = []
    rmts_super = []
    for i, j, k in zip(npts, r0s, rmts):
        for l in range(abs(det)):
            npts_super.append(i)
            r0s_super.append(j)
            rmts_super.append(k)

    w = open(filename.split('/')[-1]+"S", 'w')
    w.write( get_wien2k_struct(supercell, npts_super, r0s_super, rmts_super) )
    w.close()
    for i, cell in enumerate( cells_with_displacements ):
        symmetry = Symmetry( cell )
        print "Number of non-equivalent atoms in %sS-%03d: %d" % (
            filename, i+1, len( symmetry.get_independent_atoms() ) )
        w = open(filename.split('/')[-1]+"S-%03d" % (i+1), 'w')
        w.write( get_wien2k_struct(cell, npts_super, r0s_super, rmts_super) )
        w.close()

def gencastep(fn,modenum,basis,natom,typatsym,symprec,atpos):
    from phonopy import Phonopy
    import phonopy.structure.spglib as spg
    from phonopy.structure.atoms import PhonopyAtoms as Atoms
    from phonopy.structure.symmetry import Symmetry, find_primitive, get_pointgroup

    fh=open(fn,'w')
    unitcell = Atoms(symbols=typatsym, cell=basis, positions=atpos)
    pbasis=np.eye(3)
    for i in range(len(basis)):
        pbasis[i]=basis[i]/np.linalg.norm(basis[i])
    symmetry = Symmetry(unitcell, symprec)
    rotations = symmetry.get_symmetry_operations()['rotations']
    translations = symmetry.get_symmetry_operations()['translations']
    print('Space group International symbol: %s' % symmetry.get_international_table())

    fh.write('%BLOCK LATTICE_CART\n')
    for bl in basis:
        fh.write('%s\n' % ''.join(' %12.8f' % b for b in bl))
    fh.write('%ENDBLOCK LATTICE_CART\n\n')
    fh.write('%BLOCK POSITIONS_ABS\n')
    for i in range(len(typatsym)):
        fh.write("  %3s   " % typatsym[i])
        fh.write('%s\n' % ''.join(' %12.8f' % p for p in atpos[i].tolist()))
    fh.write('%ENDBLOCK POSITIONS_ABS\n\n')

    fh.write('SYMMETRY_TOL : %f ang\n' % symprec)
    fh.write('SYMMETRY_GENERATE \n')
    fh.write('#KPOINT_MP_GRID : 4 4 4\n#KPOINT_MP_OFFSET : 0.5 0.5 0.5\n')

def get_born_OUTCAR(
    poscar_filename="POSCAR",
    outcar_filename="OUTCAR",
    primitive_axis=np.eye(3),
    supercell_matrix=np.eye(3, dtype="intc"),
    is_symmetry=True,
    symmetrize_tensors=False,
    symprec=1e-5,
):
    ucell = read_vasp(poscar_filename)
    outcar = open(outcar_filename)

    borns, epsilon = _read_born_and_epsilon(outcar)
    num_atom = len(borns)
    assert num_atom == ucell.get_number_of_atoms()

    if symmetrize_tensors:
        lattice = ucell.get_cell().T
        positions = ucell.get_scaled_positions()
        u_sym = Symmetry(ucell, is_symmetry=is_symmetry, symprec=symprec)
        point_sym = [similarity_transformation(lattice, r) for r in u_sym.get_pointgroup_operations()]
        epsilon = _symmetrize_tensor(epsilon, point_sym)
        borns = _symmetrize_borns(borns, u_sym, lattice, positions, symprec)

    inv_smat = np.linalg.inv(supercell_matrix)
    scell = get_supercell(ucell, supercell_matrix, symprec=symprec)
    pcell = get_primitive(scell, np.dot(inv_smat, primitive_axis), symprec=symprec)
    p2s = np.array(pcell.get_primitive_to_supercell_map(), dtype="intc")
    p_sym = Symmetry(pcell, is_symmetry=is_symmetry, symprec=symprec)
    s_indep_atoms = p2s[p_sym.get_independent_atoms()]
    u2u = scell.get_unitcell_to_unitcell_map()
    u_indep_atoms = [u2u[x] for x in s_indep_atoms]
    reduced_borns = borns[u_indep_atoms].copy()

    return reduced_borns, epsilon

def _distribute_forces(supercell, disp, forces, filename, symprec):
    natom = supercell.get_number_of_atoms()
    lattice = supercell.get_cell()
    symbols = supercell.get_chemical_symbols()
    positions = supercell.get_positions()
    positions[disp[0]] += disp[1]
    cell = Atoms(cell=lattice, positions=positions, symbols=symbols, pbc=True)
    symmetry = Symmetry(cell, symprec)
    independent_atoms = symmetry.get_independent_atoms()

    # Rotation matrices in Cartesian
    rotations = []
    for r in symmetry.get_symmetry_operations()["rotations"]:
        rotations.append(similarity_transformation(lattice.T, r))

    map_operations = symmetry.get_map_operations()
    map_atoms = symmetry.get_map_atoms()

    atoms_in_dot_scf = _get_independent_atoms_in_dot_scf(filename)

    if len(forces) != len(atoms_in_dot_scf):
        print("%s does not contain necessary information." % filename)
        print('Plese check if there are "FGL" lines with')
        print('"total forces" are required.')
        return False

    if len(atoms_in_dot_scf) == natom:
        print("It is assumed that there is no symmetrically-equivalent " "atoms in ")
        print("'%s' at wien2k calculation." % filename)
        force_set = forces
    elif len(forces) != len(independent_atoms):
        print("Non-equivalent atoms of %s could not be recognized by phonopy." % filename)
        return False
    else:
        # 1. Transform wien2k forces to those on independent atoms
        indep_atoms_to_wien2k = []
        forces_remap = []
        for i, pos_wien2k in enumerate(atoms_in_dot_scf):
            for j, pos in enumerate(cell.get_scaled_positions()):
                diff = pos_wien2k - pos
                diff -= np.rint(diff)
                if (abs(diff) < symprec).all():
                    forces_remap.append(np.dot(rotations[map_operations[j]], forces[i]))
                    indep_atoms_to_wien2k.append(map_atoms[j])
                    break

        if len(forces_remap) != len(forces):
            print("Atomic position mapping between Wien2k and phonopy failed.")
            print("If you think this is caused by a bug of phonopy")
            print("please report it in the phonopy mainling list.")
            return False

        # 2. Distribute forces from independent to dependent atoms.
        force_set = []
        for i in range(natom):
            j = indep_atoms_to_wien2k.index(map_atoms[i])
            force_set.append(np.dot(rotations[map_operations[i]].T, forces_remap[j]))

    return force_set

    def test_magmom(self):
        symprec = 1e-5
        cell = get_unitcell_from_phonopy_yaml(os.path.join(data_dir,"Cr.yaml"))
        symmetry_nonspin = Symmetry(cell, symprec=symprec)
        atom_map_nonspin = symmetry_nonspin.get_map_atoms()
        len_sym_nonspin = len(
            symmetry_nonspin.get_symmetry_operations()['rotations'])
        
        spin = [1, -1]
        cell_withspin = cell.copy()
        cell_withspin.set_magnetic_moments(spin)
        symmetry_withspin = Symmetry(cell_withspin, symprec=symprec)
        atom_map_withspin = symmetry_withspin.get_map_atoms()
        len_sym_withspin = len(
            symmetry_withspin.get_symmetry_operations()['rotations'])

        broken_spin = [1, -2]
        cell_brokenspin = cell.copy()
        cell_brokenspin = cell.copy()
        cell_brokenspin.set_magnetic_moments(broken_spin)
        symmetry_brokenspin = Symmetry(cell_brokenspin, symprec=symprec)
        atom_map_brokenspin = symmetry_brokenspin.get_map_atoms()
        len_sym_brokenspin = len(
            symmetry_brokenspin.get_symmetry_operations()['rotations'])

        self.assertTrue((atom_map_nonspin == atom_map_withspin).all())
        self.assertFalse((atom_map_nonspin == atom_map_brokenspin).all())
        self.assertTrue(len_sym_nonspin == len_sym_withspin)
        self.assertFalse(len_sym_nonspin == len_sym_brokenspin)

def get_coarse_ir_grid_points(primitive,
                              mesh,
                              mesh_divisors,
                              coarse_mesh_shifts,
                              is_kappa_star=True,
                              symprec=1e-5):
    mesh = np.array(mesh, dtype='intc')

    symmetry = Symmetry(primitive, symprec)
    point_group = symmetry.get_pointgroup_operations()

    if mesh_divisors is None:
        (ir_grid_points,
         ir_grid_weights,
         grid_address,
         grid_mapping_table) = get_ir_grid_points(mesh, point_group)
    else:
        mesh_divs = np.array(mesh_divisors, dtype='intc')
        coarse_mesh = mesh // mesh_divs
        if coarse_mesh_shifts is None:
            coarse_mesh_shifts = [False, False, False]
    
        if not is_kappa_star:
            coarse_grid_address = get_grid_address(coarse_mesh)
            coarse_grid_points = np.arange(np.prod(coarse_mesh), dtype='intc')
            coarse_grid_weights = np.ones(len(coarse_grid_points), dtype='intc')
        else:
            (coarse_ir_grid_points,
             coarse_ir_grid_weights,
             coarse_grid_address,
             coarse_grid_mapping_table) = get_ir_grid_points(
                 coarse_mesh,
                 point_group,
                 mesh_shifts=coarse_mesh_shifts)
        ir_grid_points = from_coarse_to_dense_grid_points(
            mesh,
            mesh_divs,
            coarse_grid_points,
            coarse_grid_address,
            coarse_mesh_shifts=coarse_mesh_shifts)
        grid_address = get_grid_address(mesh)
        ir_grid_weights = ir_grid_weights

    primitive_lattice = np.linalg.inv(primitive.get_cell())
    bz_grid_address, bz_map = spg.relocate_BZ_grid_address(grid_address,
                                                           mesh,
                                                           primitive_lattice)

    return (ir_grid_points,
            ir_grid_weights,
            bz_grid_address,
            grid_mapping_table)

 def test_parse_wien2k_struct(self):
     cell, npts, r0s, rmts = parse_wien2k_struct("BaGa2.struct")
     lattice = cell.get_cell().T
     displacements, supercell = parse_disp_yaml("disp_BaGa2.yaml",
                                                return_cell=True)
     symmetry = Symmetry(cell)
     print(PhonopyAtoms(atoms=cell))
     sym_op = symmetry.get_symmetry_operations()
     print(symmetry.get_international_table())
     for i, (r, t) in enumerate(
             zip(sym_op['rotations'], sym_op['translations'])):
         print("--- %d ---" % (i + 1))
         print(r)
         print(t)

def get_born_parameters(f, primitive, is_symmetry):
    # Read unit conversion factor, damping factor, ...
    factors = [float(x) for x in f.readline().split()]
    if len(factors) < 1:
        print "BORN file format of line 1 is incorrect"
        return False
    if len(factors) < 2:
        factors = factors[0]

    # Read dielectric constant
    line = f.readline().split()
    if not len(line) == 9:
        print "BORN file format of line 2 is incorrect"
        return False
    dielectric = np.reshape([float(x) for x in line], (3, 3))

    # Read Born effective charge
    symmetry = Symmetry(primitive, is_symmetry=is_symmetry)
    independent_atoms = symmetry.get_independent_atoms()
    born = np.zeros((primitive.get_number_of_atoms(), 3, 3), dtype=float)

    for i in independent_atoms:
        line = f.readline().split()
        if len(line) == 0:
            print "Number of lines for Born effect charge is not enough."
            return False
        if not len(line) == 9:
            print "BORN file format of line %d is incorrect" % (i + 3)
            return False
        born[i] = np.reshape([float(x) for x in line], (3, 3))

    # Expand Born effective charges to all atoms in the primitive cell
    rotations = symmetry.get_symmetry_operations()['rotations']
    map_operations = symmetry.get_map_operations()
    map_atoms = symmetry.get_map_atoms()

    for i in range(primitive.get_number_of_atoms()):
        # R_cart = L R L^-1
        rot_cartesian = similarity_transformation(
            primitive.get_cell().transpose(), rotations[map_operations[i]])
        # R_cart^T B R_cart^-T (inverse rotation is required to transform)
        born[i] = similarity_transformation(rot_cartesian.transpose(),
                                            born[map_atoms[i]])

    non_anal = {'born': born,
                'factor': factors,
                'dielectric': dielectric }

    return non_anal

def get_number_of_triplets(primitive,
                           mesh,
                           grid_point,
                           symprec=1e-5):
    mesh = np.array(mesh, dtype='intc')
    symmetry = Symmetry(primitive, symprec)
    point_group = symmetry.get_pointgroup_operations()
    primitive_lattice = np.linalg.inv(primitive.get_cell())
    triplets_at_q, _, _, _, _, _ = get_triplets_at_q(
        grid_point,
        mesh,
        point_group,
        primitive_lattice)

    return len(triplets_at_q)

def parse_BORN( primitive,  filename = "BORN" ):
    file = open( filename, 'r' )

    # Read unit conversion factor, damping factor, ...
    factors = [ float( x ) for x in file.readline().split() ]
    if len( factors ) < 1:
        print "BORN file format of line 1 is incorrect"
        return None
    if len( factors ) < 2:
        factors.append( Damping_Factor )

    # Read dielectric constant
    line = file.readline().split()
    if 9 < len( line ) or len( line ) < 9:
        print "BORN file format of line 2 is incorrect"
        return None
    dielectric = np.reshape( [ float( x ) for x in line ], ( 3, 3 ) )

    # Read Born effective charge
    symmetry = Symmetry( primitive )
    independent_atoms = symmetry.get_independent_atoms()
    born = np.zeros( ( primitive.get_number_of_atoms(), 3, 3 ), dtype=float )

    for i in independent_atoms:
        line = file.readline().split()
        if 9 < len( line ) or len( line ) < 9:
            print "BORN file format of line %d is incorrect" % ( i + 3 )
            return None
        born[ i ] = np.reshape( [ float( x ) for x in line ], ( 3, 3 ) )

    # Expand Born effective charges to all atoms in the primitive cell
    rotations = symmetry.get_symmetry_operations()['rotations']
    map_operations = symmetry.get_map_operations()
    map_atoms = symmetry.get_map_atoms()

    for i in range( primitive.get_number_of_atoms() ):
        # R_cart = L R L^-1
        rot_cartesian = similarity_transformation(
            primitive.get_cell().transpose(), rotations[ map_operations[i] ] )
        # R_cart^T B R_cart^-T ( inverse rotation is required to transform )
        born[i] = similarity_transformation( rot_cartesian.transpose(),
                                             born[ map_atoms[ i ] ] )

    non_anal = {'born': born,
                'factor': factors,
                'dielectric': dielectric }

    return non_anal

 def _search_primitive_symmetry(self):
     self._primitive_symmetry = Symmetry(self._primitive,
                                         self._symprec,
                                         self._is_symmetry)
     if (len(self._symmetry.get_pointgroup_operations()) !=
         len(self._primitive_symmetry.get_pointgroup_operations())):
         print("Warning: point group symmetries of supercell and primitive"
               "cell are different.")

    def __init__(self,
                 mesh,
                 primitive,
                 supercell,
                 fc2,
                 nac_params=None,
                 nac_q_direction=None,
                 sigma=None,
                 cutoff_frequency=None,
                 frequency_step=None,
                 num_frequency_points=None,
                 temperatures=None,
                 frequency_factor_to_THz=VaspToTHz,
                 frequency_scale_factor=1.0,
                 is_nosym=False,
                 symprec=1e-5,
                 filename=None,
                 log_level=False,
                 lapack_zheev_uplo='L'):

        self._grid_point = None
        self._mesh = np.array(mesh, dtype='intc')
        self._primitive = primitive
        self._supercell = supercell
        self._fc2 = fc2
        self._nac_params = nac_params
        self._nac_q_direction = None
        self.set_nac_q_direction(nac_q_direction)
        self._sigma = None
        self.set_sigma(sigma)

        if cutoff_frequency is None:
            self._cutoff_frequency = 0
        else:
            self._cutoff_frequency = cutoff_frequency
        self._frequency_step = frequency_step
        self._num_frequency_points = num_frequency_points
        self._temperatures = temperatures
        self._frequency_factor_to_THz = frequency_factor_to_THz
        self._frequency_scale_factor = frequency_scale_factor
        self._is_nosym = is_nosym
        self._symprec = symprec
        self._filename = filename
        self._log_level = log_level
        self._lapack_zheev_uplo = lapack_zheev_uplo

        self._num_band = self._primitive.get_number_of_atoms() * 3
        self._reciprocal_lattice = np.linalg.inv(self._primitive.get_cell())
        self._set_dynamical_matrix()
        self._symmetry = Symmetry(primitive, symprec)

        self._tetrahedron_method = None
        self._phonon_done = None
        self._frequencies = None
        self._eigenvectors = None
            
        self._joint_dos = None
        self._frequency_points = None