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

Abinit support

This commit is contained in:
Atsushi Togo 2014-10-29 17:18:56 +09:00
parent cdb171ebf1
commit 86c96c90e9
4 changed files with 111 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
phonopy/file_IO.py
View File

@ -36,44 +36,38 @@ import sys
import StringIO
import numpy as np
import phonopy.interface.vasp as vasp
import phonopy.interface.wien2k as wien2k
from phonopy.structure.symmetry import Symmetry
from phonopy.harmonic.force_constants import similarity_transformation
from phonopy.structure.atoms import Atoms
from phonopy.cui.settings import fracval
# Constant
Damping_Factor = 0.25
# Utility to read file with ignoring blank lines
def get_line_ignore_blank(f):
def _get_line_ignore_blank(f):
line = f.readline().strip()
if line == '':
line = get_line_ignore_blank(f)
line = _get_line_ignore_blank(f)
return line
# Parse FORCE_SETS
def parse_FORCE_SETS(is_translational_invariance=False, filename="FORCE_SETS"):
f = open(filename, 'r')
return get_set_of_forces(f, is_translational_invariance)
return _get_set_of_forces(f, is_translational_invariance)
def parse_FORCE_SETS_from_strings(strings, is_translational_invariance=False):
return get_set_of_forces(StringIO.StringIO(strings),
return _get_set_of_forces(StringIO.StringIO(strings),
is_translational_invariance)
def get_set_of_forces(f, is_translational_invariance):
def _get_set_of_forces(f, is_translational_invariance):
set_of_forces = []
num_atom = int(get_line_ignore_blank(f))
num_displacements = int(get_line_ignore_blank(f))
num_atom = int(_get_line_ignore_blank(f))
num_displacements = int(_get_line_ignore_blank(f))
for i in range(num_displacements):
line = get_line_ignore_blank(f)
line = _get_line_ignore_blank(f)
atom_number = int(line)
line = get_line_ignore_blank(f).split()
line = _get_line_ignore_blank(f).split()
displacement = np.array([float(x) for x in line])
forces_tmp = []
for j in range(num_atom):
line = get_line_ignore_blank(f).split()
line = _get_line_ignore_blank(f).split()
forces_tmp.append(np.array([float(x) for x in line]))
forces_tmp = np.array(forces_tmp, dtype='double')
@ -90,7 +84,6 @@ def get_set_of_forces(f, is_translational_invariance):
return dataset
# Parse QPOINTS
def parse_QPOINTS(filename="QPOINTS"):
f = open(filename, 'r')
num_qpoints = int(f.readline().strip())
@ -99,14 +92,42 @@ def parse_QPOINTS(filename="QPOINTS"):
qpoints.append([fracval(x) for x in f.readline().strip().split()])
return np.array(qpoints)
# Write FORCE_SETS for VASP
def _get_drift_forces(forces, filename=None):
drift_force = np.sum(forces, axis=0) / len(forces)
if filename is None:
print "Drift force"
else:
print "Drift force of %s" % filename
print "%12.8f %12.8f %12.8f" % tuple(drift_force)
print "This drift force was subtracted from forces."
return drift_force
def write_FORCE_SETS_abinit(forces_filenames,
displacements,
num_atom,
filename='FORCE_SETS'):
import phonopy.interface.abinit as abinit
for abinit_filename, disp in zip(forces_filenames,
displacements['first_atoms']):
abinit_forces = abinit.get_forces_abinit(abinit_filename, num_atom)
if abinit_forces is False:
return False
drift_force = _get_drift_forces(abinit_forces)
disp['forces'] = np.array(abinit_forces) - drift_force
write_FORCE_SETS(displacements, filename=filename)
return True
def write_FORCE_SETS_wien2k(forces_filenames,
displacements,
supercell,
filename='FORCE_SETS',
is_zero_point=False,
is_distribute=True,
symprec=1e-5):
import phonopy.interface.wien2k as wien2k
natom = supercell.get_number_of_atoms()
lattice = supercell.get_cell()
@ -132,12 +153,7 @@ def write_FORCE_SETS_wien2k(forces_filenames,
else:
force_set = wien2k_forces
drift_force = np.sum(force_set, axis=0) / len(force_set)
print "Drift force of %s" % wien2k_filename
print "%12.8f %12.8f %12.8f" % tuple(drift_force)
print "This drift force was subtracted from forces."
print
drift_force = _get_drift_forces(force_set, filename=wien2k_filename)
disp['forces'] = np.array(force_set) - drift_force
write_FORCE_SETS(displacements, filename=filename)
@ -182,7 +198,7 @@ def write_FORCE_SETS_vasp(forces_filenames,
print "%d" % (count + 1),
count += 1
if not check_forces(zero_forces, num_atom, forces_filenames[0]):
if not _check_forces(zero_forces, num_atom, forces_filenames[0]):
are_files_correct = False
else:
force_files = forces_filenames
@ -200,7 +216,7 @@ def write_FORCE_SETS_vasp(forces_filenames,
print "%d" % (count + 1),
count += 1
if not check_forces(disp['forces'], num_atom, force_files[i]):
if not _check_forces(disp['forces'], num_atom, force_files[i]):
are_files_correct = False
if verbose:
@ -212,7 +228,7 @@ def write_FORCE_SETS_vasp(forces_filenames,
return are_files_correct
def check_forces(forces, num_atom, filename):
def _check_forces(forces, num_atom, filename):
if len(forces) != num_atom:
print " \"%s\" does not contain necessary information." % filename,
return False
@ -239,9 +255,6 @@ def write_FORCE_SETS(dataset, filename='FORCE_SETS', zero_forces=None):
for f in forces[count]:
fp.write("%15.10f %15.10f %15.10f\n" % (tuple(f)))
def mycmp(x, y):
return cmp(x[0], y[0])
def parse_disp_yaml(filename="disp.yaml", return_cell=False):
try:
import yaml

25
phonopy/interface/abinit.py
View File

@ -58,6 +58,31 @@ def write_abinit(filename, cell):
f = open(filename, 'w')
f.write(get_abinit_structure(cell))
def write_supercells_with_displacements(supercell,
cells_with_displacements):
write_abinit("supercell.in", supercell)
for i, cell in enumerate(cells_with_displacements):
write_abinit("supercell-%03d.in" % (i + 1), cell)
def get_forces_abinit(filename, num_atom):
f = open(filename)
for line in f:
if 'cartesian forces (eV/Angstrom)' in line:
break
forces = []
for line in f:
elems = line.split()
if len(elems) > 3:
forces.append([float(x) for x in elems[1:4]])
else:
return False
if len(forces) == num_atom:
break
return forces
def get_abinit_structure(cell):
znucl = []
numbers = cell.get_atomic_numbers()

1
phonopy/units.py
View File

@ -62,4 +62,5 @@ VaspToTHz = sqrt(EV/AMU)/Angstrom/(2*pi)/1e12 # [THz] 15.633302
VaspToCm = VaspToTHz * THzToCm # [cm^-1] 521.47083
EvTokJmol = EV / 1000 * Avogadro # [kJ/mol] 96.4853910
Wien2kToTHz = sqrt(Rydberg/1000*EV/AMU)/(Bohr*1e-10)/(2*pi)/1e12 # [THz] 3.44595837
AbinitToTHz = sqrt(EV/(AMU*Bohr))/Angstrom/(2*pi)/1e12 # [THz] 15.633302
EVAngstromToGPa = EV * 1e21

53
scripts/phonopy
View File

@ -40,8 +40,6 @@ import numpy as np
from phonopy import *
from optparse import OptionParser
import phonopy.file_IO as file_IO
import phonopy.interface.wien2k as wien2k
import phonopy.interface.vasp as vasp
from phonopy.cui.settings import Settings, PhonopyConfParser
from phonopy.cui.show_symmetry import check_symmetry
from phonopy.units import *
@ -104,15 +102,15 @@ def file_exists(filename, log_level):
parser = OptionParser()
parser.set_defaults(
displacement_distance=None,
abinit_struct_file=None,
anime=None,
band_indices=None,
band_paths=None,
band_points=None,
cutoff_frequency=None,
cutoff_radius=None,
displacement_distance=None,
dynamical_matrix_decimals=None,
is_dos_mode=False,
factor=None,
fc_symmetry=None,
fc_computation_algorithm=None,
@ -126,6 +124,7 @@ parser.set_defaults(
is_band_connection=False,
is_check_symmetry=False,
is_displacement=False,
is_dos_mode=False,
is_eigenvectors=False,
is_gamma_center=False,
is_graph_plot=False,
@ -168,6 +167,9 @@ parser.set_defaults(
verbose=False,
wien2k_struct_file=None)
parser.add_option("--abinit", dest="abinit_struct_file",
action="store", type="string",
help="Read Abinit structure from .in file", metavar="FILE")
parser.add_option("--amplitude", dest="displacement_distance", type="float",
help="Distance of displacements")
parser.add_option("--anime", dest="anime",
@ -378,8 +380,13 @@ if log_level > 0:
# Phonopy interface mode
if options.wien2k_struct_file:
interface_mode = 'wien2k'
import phonopy.interface.wien2k as wien2k
elif options.abinit_struct_file:
interface_mode = 'abinit'
import phonopy.interface.abinit as abinit
else:
interface_mode = 'vasp'
import phonopy.interface.vasp as vasp
######################################################
# Create FORCE_CONSTANTS (--fc or --force_constants) #
@ -424,7 +431,7 @@ if options.force_sets_mode or options.force_sets_zero_mode:
if file_exists('disp.yaml', log_level):
if interface_mode == 'vasp':
displacements = file_IO.parse_disp_yaml('disp.yaml')
if interface_mode == 'wien2k':
if interface_mode == 'wien2k' or interface_mode == 'abinit':
displacements, supercell = file_IO.parse_disp_yaml(
'disp.yaml', return_cell=True)
@ -445,9 +452,21 @@ if options.force_sets_mode or options.force_sets_zero_mode:
is_created = file_IO.write_FORCE_SETS_vasp(
args,
displacements,
'FORCE_SETS',
filename='FORCE_SETS',
is_zero_point=options.force_sets_zero_mode)
if interface_mode == 'abinit':
print "**********************************************************"
print "**** Abinit FORCE_SETS support is experimental. ****"
print "**** Your feedback would be appreciated. ****"
print "**********************************************************"
is_created = file_IO.write_FORCE_SETS_abinit(
args,
displacements,
supercell.get_number_of_atoms(),
filename='FORCE_SETS')
if interface_mode == 'wien2k':
print "**********************************************************"
print "**** Wien2k FORCE_SETS support is experimental. ****"
@ -457,8 +476,7 @@ if options.force_sets_mode or options.force_sets_zero_mode:
args,
displacements,
supercell,
'FORCE_SETS',
is_zero_point=options.force_sets_zero_mode,
filename='FORCE_SETS',
is_distribute=(not options.is_wien2k_p1),
symprec=options.symprec)
@ -473,7 +491,6 @@ if options.force_sets_mode or options.force_sets_zero_mode:
##########################
# Read crystal structure #
##########################
# Parse VASP-type crystal structure
if interface_mode == 'vasp':
if options.cell_poscar is None:
if file_exists("POSCAR", log_level):
@ -484,7 +501,10 @@ if interface_mode == 'vasp':
unitcell = vasp.read_vasp(unitcell_filename)
# Parse Wien2k-type crystal structure
if interface_mode == 'abinit':
if file_exists(options.abinit_struct_file, log_level):
unitcell = abinit.read_abinit(options.abinit_struct_file)
if interface_mode == 'wien2k':
if file_exists(options.wien2k_struct_file, log_level):
unitcell, npts, r0s, rmts = wien2k.parse_wien2k_struct(
@ -529,13 +549,16 @@ else:
run_mode = settings.get_run_mode()
# Physical units: energy, distance, atomic mass
# vasp : eV, Angstrom, AMU
# wien2k : Ry, au, AMU
# Physical units: energy, distance, atomic mass, force
# vasp : eV, Angstrom, AMU, eV/Angstrom
# wien2k : Ry, au, AMU, Ry/au
# abinit : hartree, au(=bohr), AMU, eV/Angstrom
# --factor: Frequency conversion factor
if options.factor is None:
if interface_mode == 'wien2k':
factor = Wien2kToTHz
elif interface_mode == 'abinit':
factor = AbinitToTHz
else:
factor = VaspToTHz
else:
@ -543,7 +566,7 @@ else:
# --amplitude
if settings.get_displacement_distance() is None:
if interface_mode == 'wien2k':
if interface_mode == 'wien2k' or interface_mode == 'abinit':
displacement_distance = 0.02
else:
displacement_distance = 0.01
@ -775,6 +798,8 @@ if run_mode == 'displacements':
rmts,
settings.get_supercell_matrix(),
filename=options.wien2k_struct_file)
elif interface_mode == 'abinit':
abinit.write_supercells_with_displacements(supercell, cells_with_disps)
else:
vasp.write_supercells_with_displacements(supercell, cells_with_disps)