mirror of https://github.com/phonopy/phonopy.git
Refactoring and some test for wien2k interface
This commit is contained in:
Atsushi Togo
parent
3256741b6a
commit
f0078539d5
|
|
@ -52,7 +52,7 @@ def parse_set_of_forces(displacements,
|
|||
# Parse wien2k case.scf file
|
||||
wien2k_forces = get_forces_wien2k(wien2k_filename, lattice)
|
||||
if is_distribute:
|
||||
force_set = distribute_forces(
|
||||
force_set = _distribute_forces(
|
||||
supercell,
|
||||
[disp['number'], disp['displacement']],
|
||||
wien2k_forces,
|
||||
|
|
@ -62,8 +62,8 @@ def parse_set_of_forces(displacements,
|
|||
return False
|
||||
else:
|
||||
if not (natom == len(wien2k_forces)):
|
||||
print "%s contains only forces of %d atoms" % (
|
||||
wien2k_filename, len(wien2k_forces))
|
||||
print("%s contains only forces of %d atoms" %
|
||||
(wien2k_filename, len(wien2k_forces)))
|
||||
return False
|
||||
else:
|
||||
force_set = wien2k_forces
|
||||
|
|
@ -87,7 +87,7 @@ def create_FORCE_SETS(forces_filenames,
|
|||
# Parse wien2k case.scf file
|
||||
wien2k_forces = get_forces_wien2k(wien2k_filename, lattice)
|
||||
if is_distribute:
|
||||
force_set = distribute_forces(
|
||||
force_set = _distribute_forces(
|
||||
supercell,
|
||||
[disp['number'], disp['displacement']],
|
||||
wien2k_forces,
|
||||
|
|
@ -97,8 +97,8 @@ def create_FORCE_SETS(forces_filenames,
|
|||
return False
|
||||
else:
|
||||
if not (natom == len(wien2k_forces)):
|
||||
print "%s contains only forces of %d atoms" % (
|
||||
wien2k_filename, len(wien2k_forces))
|
||||
print("%s contains only forces of %d atoms" %
|
||||
(wien2k_filename, len(wien2k_forces)))
|
||||
return False
|
||||
else:
|
||||
force_set = wien2k_forces
|
||||
|
|
@ -298,9 +298,10 @@ def write_supercells_with_displacements(supercell,
|
|||
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')
|
||||
supercell_filename = filename.split('/')[-1]+"S-%03d" % (i + 1)
|
||||
print("Number of non-equivalent atoms in %s: %d" %
|
||||
(supercell_filename, len(symmetry.get_independent_atoms())))
|
||||
w = open(supercell_filename, 'w')
|
||||
w.write(get_wien2k_struct(cell, npts_super, r0s_super, rmts_super))
|
||||
w.close()
|
||||
|
||||
|
|
@ -323,24 +324,7 @@ def get_forces_wien2k(filename, lattice):
|
|||
|
||||
return forces[-num_atom:]
|
||||
|
||||
def get_independent_atoms_in_dot_scf(filename):
|
||||
positions = []
|
||||
for line in open(filename):
|
||||
if line[:4] == ":POS":
|
||||
if "POSITION" in line:
|
||||
x = float(line[30:37])
|
||||
y = float(line[38:45])
|
||||
z = float(line[46:53])
|
||||
else:
|
||||
x = float(line[27:34])
|
||||
y = float(line[35:42])
|
||||
z = float(line[43:50])
|
||||
num_atom = int(line[4:7])
|
||||
positions.append([x,y,z])
|
||||
|
||||
return np.array(positions)[-num_atom:]
|
||||
|
||||
def distribute_forces(supercell, disp, forces, filename, symprec):
|
||||
def _distribute_forces(supercell, disp, forces, filename, symprec):
|
||||
natom = supercell.get_number_of_atoms()
|
||||
lattice = supercell.get_cell()
|
||||
symbols = supercell.get_chemical_symbols()
|
||||
|
|
@ -361,21 +345,24 @@ def distribute_forces(supercell, disp, forces, filename, symprec):
|
|||
map_operations = symmetry.get_map_operations()
|
||||
map_atoms = symmetry.get_map_atoms()
|
||||
|
||||
atoms_in_dot_scf = get_independent_atoms_in_dot_scf(filename)
|
||||
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."
|
||||
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
|
||||
print ""
|
||||
print('')
|
||||
print("It is assumed that there is no symmetrically-equivalent "
|
||||
"atoms in ")
|
||||
print("\'%s\' at wien2k calculation." % filename)
|
||||
print('')
|
||||
force_set = forces
|
||||
elif len(forces) != len(independent_atoms):
|
||||
print "Non-equivalent atoms of %s could not be recognized by phonopy." % filename
|
||||
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
|
||||
|
|
@ -385,27 +372,46 @@ def distribute_forces(supercell, disp, forces, filename, symprec):
|
|||
for j, pos in enumerate(cell.get_scaled_positions()):
|
||||
diff = pos_wien2k - pos
|
||||
diff -= np.rint(diff)
|
||||
if (abs(diff) < 0.00001).all():
|
||||
if (abs(diff) < symprec).all():
|
||||
forces_remap.append(
|
||||
np.dot(forces[i], rotations[map_operations[j]].T))
|
||||
np.dot(rotations[map_operations[j]], forces[i]))
|
||||
indep_atoms_to_wien2k.append(map_atoms[j])
|
||||
break
|
||||
|
||||
if not 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."
|
||||
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(
|
||||
forces_remap[indep_atoms_to_wien2k.index(map_atoms[i])],
|
||||
rotations[map_operations[i]]))
|
||||
rotations[map_operations[i]].T, forces_remap[j]))
|
||||
|
||||
return force_set
|
||||
|
||||
def _get_independent_atoms_in_dot_scf(filename):
|
||||
positions = []
|
||||
for line in open(filename):
|
||||
if line[:4] == ":POS":
|
||||
if "POSITION" in line:
|
||||
x = float(line[30:37])
|
||||
y = float(line[38:45])
|
||||
z = float(line[46:53])
|
||||
else:
|
||||
x = float(line[27:34])
|
||||
y = float(line[35:42])
|
||||
z = float(line[43:50])
|
||||
num_atom = int(line[4:7])
|
||||
positions.append([x,y,z])
|
||||
|
||||
return np.array(positions)[-num_atom:]
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from optparse import OptionParser
|
||||
from phonopy.interface.vasp import write_vasp, read_vasp
|
||||
|
|
@ -437,7 +443,7 @@ if __name__ == '__main__':
|
|||
cell.set_cell(lattice)
|
||||
npts, r0s, rmts = parse_core_param(open(args[1]))
|
||||
text = get_wien2k_struct(cell, npts, r0s, rmts)
|
||||
print text
|
||||
print(text)
|
||||
|
||||
elif options.w2v:
|
||||
cell, npts, r0s, rmts = parse_wien2k_struct(args[0])
|
||||
|
|
@ -455,4 +461,5 @@ if __name__ == '__main__':
|
|||
w.write("%-10s %5d %10.8f %10.5f\n" %
|
||||
(symbol, npt, r0, rmt))
|
||||
else:
|
||||
print "You need to set -r or -w option."
|
||||
print("You need to set -r or -w option.")
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,114 @@
|
|||
blebleble
|
||||
P LATTICE,NONEQUIV.ATOMS: 2 191_P6/mmm
|
||||
MODE OF CALC=RELA unit=bohr
|
||||
8.429120 8.429120 9.701070 90.000000 90.000000120.000000
|
||||
ATOM -1: X=0.00000000 Y=0.00000000 Z=0.00000000
|
||||
MULT= 1 ISPLIT= 4
|
||||
Ba0+ NPT= 781 R0=0.00001000 RMT= 2.5000 Z: 56.0
|
||||
LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000
|
||||
0.0000000 1.0000000 0.0000000
|
||||
0.0000000 0.0000000 1.0000000
|
||||
ATOM -2: X=0.33333333 Y=0.66666667 Z=0.50000000
|
||||
MULT= 2 ISPLIT= 4
|
||||
-2: X=0.66666667 Y=0.33333333 Z=0.50000000
|
||||
Ga0+ NPT= 781 R0=0.00005000 RMT= 2.2800 Z: 31.0
|
||||
LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000
|
||||
0.0000000 1.0000000 0.0000000
|
||||
0.0000000 0.0000000 1.0000000
|
||||
24 NUMBER OF SYMMETRY OPERATIONS
|
||||
-1 0 0 0.00000000
|
||||
-1 1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
1
|
||||
-1 0 0 0.00000000
|
||||
-1 1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
2
|
||||
-1 1 0 0.00000000
|
||||
-1 0 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
3
|
||||
-1 1 0 0.00000000
|
||||
-1 0 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
4
|
||||
-1 0 0 0.00000000
|
||||
0-1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
5
|
||||
-1 0 0 0.00000000
|
||||
0-1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
6
|
||||
-1 1 0 0.00000000
|
||||
0 1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
7
|
||||
-1 1 0 0.00000000
|
||||
0 1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
8
|
||||
0-1 0 0.00000000
|
||||
-1 0 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
9
|
||||
0-1 0 0.00000000
|
||||
-1 0 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
10
|
||||
0 1 0 0.00000000
|
||||
-1 1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
11
|
||||
0 1 0 0.00000000
|
||||
-1 1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
12
|
||||
0-1 0 0.00000000
|
||||
1-1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
13
|
||||
0-1 0 0.00000000
|
||||
1-1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
14
|
||||
0 1 0 0.00000000
|
||||
1 0 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
15
|
||||
0 1 0 0.00000000
|
||||
1 0 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
16
|
||||
1-1 0 0.00000000
|
||||
0-1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
17
|
||||
1-1 0 0.00000000
|
||||
0-1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
18
|
||||
1 0 0 0.00000000
|
||||
0 1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
19
|
||||
1 0 0 0.00000000
|
||||
0 1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
20
|
||||
1-1 0 0.00000000
|
||||
1 0 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
21
|
||||
1-1 0 0.00000000
|
||||
1 0 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
22
|
||||
1 0 0 0.00000000
|
||||
1-1 0 0.00000000
|
||||
0 0-1 0.00000000
|
||||
23
|
||||
1 0 0 0.00000000
|
||||
1-1 0 0.00000000
|
||||
0 0 1 0.00000000
|
||||
24
|
||||
|
|
@ -0,0 +1,65 @@
|
|||
natom: 24
|
||||
displacements:
|
||||
- atom: 1
|
||||
direction:
|
||||
[ 1.0000000000000000, 0.0000000000000000, 1.0000000000000000 ]
|
||||
displacement:
|
||||
[ 0.0227205761666499, -0.0131177307659587, 0.0301943795797708 ]
|
||||
- atom: 9
|
||||
direction:
|
||||
[ 1.0000000000000000, 0.0000000000000000, 1.0000000000000000 ]
|
||||
displacement:
|
||||
[ 0.0227205761666499, -0.0131177307659587, 0.0301943795797708 ]
|
||||
lattice:
|
||||
- [ 14.599664103094977, -8.429119999999996, 0.000000000000001 ]
|
||||
- [ 0.000000000000000, 16.858239999999999, 0.000000000000001 ]
|
||||
- [ 0.000000000000000, 0.000000000000000, 19.402139999999999 ]
|
||||
atoms:
|
||||
- symbol: Ba # 1
|
||||
position: [ 0.00000000000000, 0.00000000000000, 0.00000000000000 ]
|
||||
- symbol: Ba # 2
|
||||
position: [ 0.50000000000000, 0.00000000000000, 0.00000000000000 ]
|
||||
- symbol: Ba # 3
|
||||
position: [ 0.00000000000000, 0.50000000000000, 0.00000000000000 ]
|
||||
- symbol: Ba # 4
|
||||
position: [ 0.50000000000000, 0.50000000000000, 0.00000000000000 ]
|
||||
- symbol: Ba # 5
|
||||
position: [ 0.00000000000000, 0.00000000000000, 0.50000000000000 ]
|
||||
- symbol: Ba # 6
|
||||
position: [ 0.50000000000000, 0.00000000000000, 0.50000000000000 ]
|
||||
- symbol: Ba # 7
|
||||
position: [ 0.00000000000000, 0.50000000000000, 0.50000000000000 ]
|
||||
- symbol: Ba # 8
|
||||
position: [ 0.50000000000000, 0.50000000000000, 0.50000000000000 ]
|
||||
- symbol: Ga # 9
|
||||
position: [ 0.16666666500000, 0.33333333500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 10
|
||||
position: [ 0.66666666500000, 0.33333333500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 11
|
||||
position: [ 0.16666666500000, 0.83333333500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 12
|
||||
position: [ 0.66666666500000, 0.83333333500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 13
|
||||
position: [ 0.16666666500000, 0.33333333500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 14
|
||||
position: [ 0.66666666500000, 0.33333333500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 15
|
||||
position: [ 0.16666666500000, 0.83333333500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 16
|
||||
position: [ 0.66666666500000, 0.83333333500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 17
|
||||
position: [ 0.33333333500000, 0.16666666500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 18
|
||||
position: [ 0.83333333500000, 0.16666666500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 19
|
||||
position: [ 0.33333333500000, 0.66666666500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 20
|
||||
position: [ 0.83333333500000, 0.66666666500000, 0.25000000000000 ]
|
||||
- symbol: Ga # 21
|
||||
position: [ 0.33333333500000, 0.16666666500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 22
|
||||
position: [ 0.83333333500000, 0.16666666500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 23
|
||||
position: [ 0.33333333500000, 0.66666666500000, 0.75000000000000 ]
|
||||
- symbol: Ga # 24
|
||||
position: [ 0.83333333500000, 0.66666666500000, 0.75000000000000 ]
|
||||
|
|
@ -0,0 +1,35 @@
|
|||
import unittest
|
||||
|
||||
import numpy as np
|
||||
from phonopy.structure.atoms import PhonopyAtoms
|
||||
from phonopy.structure.symmetry import Symmetry
|
||||
from phonopy.interface.wien2k import parse_wien2k_struct
|
||||
from phonopy.file_IO import parse_disp_yaml
|
||||
|
||||
class TestWien2k(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
|
||||
def tearDown(self):
|
||||
pass
|
||||
|
||||
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)
|
||||
|
||||
if __name__ == '__main__':
|
||||
suite = unittest.TestLoader().loadTestsFromTestCase(TestWien2k)
|
||||
unittest.TextTestRunner(verbosity=2).run(suite)
|
||||
# unittest.main()
|
||||
Loading…
Reference in New Issue