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

write_FORCE_SETS is modified to recieve displacement_dataset

This commit is contained in:
Atsushi Togo 2014-02-16 17:54:48 +09:00
parent 3fc3bbc6b0
commit dcce5fe198
2 changed files with 104 additions and 135 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

224
phonopy/file_IO.py
View File

@ -109,25 +109,26 @@ def write_FORCE_SETS_wien2k(forces_filenames,
is_distribute=True, is_distribute=True,
symprec=1e-5): symprec=1e-5):
forces = []
natom = supercell.get_number_of_atoms() natom = supercell.get_number_of_atoms()
lattice = supercell.get_cell() lattice = supercell.get_cell()
for wien2k_filename, disp in zip(forces_filenames, displacements): for wien2k_filename, disp in zip(forces_filenames,
displacements['first_atoms']):
# Parse wien2k case.scf file # Parse wien2k case.scf file
wien2k_forces = wien2k.get_forces_wien2k(wien2k_filename, lattice) wien2k_forces = wien2k.get_forces_wien2k(wien2k_filename, lattice)
if is_distribute: if is_distribute:
force_set = wien2k.distribute_forces(supercell, force_set = wien2k.distribute_forces(
disp, supercell,
wien2k_forces, [disp['number'], disp['displacement']],
wien2k_filename, wien2k_forces,
symprec) wien2k_filename,
symprec)
if not force_set: if not force_set:
return False return False
else: else:
if not (natom == len(wien2k_forces)): if not (natom == len(wien2k_forces)):
print "%s contains only forces of %d atoms" % (wien2k_filename, print "%s contains only forces of %d atoms" % (
len(wien2k_forces)) wien2k_filename, len(wien2k_forces))
return False return False
else: else:
force_set = wien2k_forces force_set = wien2k_forces
@ -136,26 +137,31 @@ def write_FORCE_SETS_wien2k(forces_filenames,
print "Drift force of %s" % wien2k_filename print "Drift force of %s" % wien2k_filename
print "%12.8f %12.8f %12.8f" % tuple(drift_force) print "%12.8f %12.8f %12.8f" % tuple(drift_force)
print "This drift force was subtracted from forces." print "This drift force was subtracted from forces."
print print
forces.append(np.array(force_set) - drift_force) disp['forces'] = np.array(force_set) - drift_force
write_FORCE_SETS(filename, natom, displacements, forces) write_FORCE_SETS(displacements, filename=filename)
return True return True
def write_FORCE_SETS_vasp(forces_filenames, def write_FORCE_SETS_vasp(forces_filenames,
displacements, displacements,
num_atom,
filename='FORCE_SETS', filename='FORCE_SETS',
is_zero_point=False, is_zero_point=False,
verbose=True): verbose=True):
try: try:
from lxml import etree from lxml import etree
except ImportError: except ImportError:
print "You need to install python-lxml." print "You need to install python-lxml."
sys.exit(1) sys.exit(1)
if verbose:
print "counter (file index):",
num_atom = displacements['natom']
count = 0
are_files_correct = True
if is_zero_point: if is_zero_point:
force_files = forces_filenames[1:] force_files = forces_filenames[1:]
if vasp.is_version528(forces_filenames[0]): if vasp.is_version528(forces_filenames[0]):
@ -164,103 +170,72 @@ def write_FORCE_SETS_vasp(forces_filenames,
else: else:
zero_forces = vasp.get_forces_vasprun_xml( zero_forces = vasp.get_forces_vasprun_xml(
etree.iterparse(forces_filenames[0], tag='varray')) etree.iterparse(forces_filenames[0], tag='varray'))
if verbose:
print "%d" % (count + 1),
count += 1
if not check_forces(zero_forces, num_atom, forces_filenames[0]):
are_files_correct = False
else: else:
force_files = forces_filenames force_files = forces_filenames
zero_forces = None zero_forces = None
forces = [] for i, disp in enumerate(displacements['first_atoms']):
# Show progress
for i in range(len(displacements)):
if vasp.is_version528(force_files[i]): if vasp.is_version528(force_files[i]):
forces.append(vasp.get_forces_vasprun_xml(etree.iterparse( disp['forces'] = vasp.get_forces_vasprun_xml(etree.iterparse(
vasp.VasprunWrapper(force_files[i]), tag='varray'))) vasp.VasprunWrapper(force_files[i]), tag='varray'))
else: else:
forces.append(vasp.get_forces_vasprun_xml( disp['forces'] = vasp.get_forces_vasprun_xml(
etree.iterparse(force_files[i], tag='varray'))) etree.iterparse(force_files[i], tag='varray'))
if is_zero_point:
dummy_forces = [zero_forces] + forces
else:
dummy_forces = forces
if is_forces_read(dummy_forces, num_atom, forces_filenames):
if verbose: if verbose:
print >> sys.stderr, "counter (file index):", print "%d" % (count + 1),
write_FORCE_SETS(filename, count += 1
num_atom,
displacements, if not check_forces(disp['forces'], num_atom, force_files[i]):
forces, are_files_correct = False
zero_forces,
verbose) if verbose:
return True print
else:
write_FORCE_SETS(displacements,
filename=filename,
zero_forces=zero_forces)
return are_files_correct
def check_forces(forces, num_atom, filename):
if len(forces) != num_atom:
print " \"%s\" does not contain necessary information." % filename,
return False return False
else:
return True
def is_forces_read(force_sets, num_atom, filenames): def write_FORCE_SETS(dataset, filename='FORCE_SETS', zero_forces=None):
is_read = True num_atom = dataset['natom']
for i, forces in enumerate(force_sets): displacements = dataset['first_atoms']
if not len(forces)==num_atom:
is_read = False
print "\'%s\' does not contain necessary information." % filenames[i]
return is_read
def write_FORCE_SETS_from_dataset(dataset,
filename='FORCE_SETS',
verbose=False):
displacements = [[x['number'], x['displacement']]
for x in dataset['first_atoms']]
forces = [x['forces'] for x in dataset['first_atoms']] forces = [x['forces'] for x in dataset['first_atoms']]
write_FORCE_SETS(filename,
dataset['natom'],
displacements,
forces,
verbose=verbose)
def write_FORCE_SETS(filename,
natom,
displacements,
forces,
zero_forces=None,
verbose=True):
disps = sort_displacements(displacements)
# Write FORCE_SETS # Write FORCE_SETS
file = open(filename, 'w') fp = open(filename, 'w')
file.write("%-5d\n" % natom) fp.write("%-5d\n" % num_atom)
file.write("%-5d\n" % len(disps)) fp.write("%-5d\n" % len(displacements))
for count, disp in enumerate(disps): for count, disp in enumerate(displacements):
# Show progress fp.write("\n%-5d\n" % (disp['number'] + 1))
if verbose: fp.write("%20.16f %20.16f %20.16f\n" % (tuple(disp['displacement'])))
print >> sys.stderr, "%d (%d) " % (count+1, disp[2] + 1),
file.write("\n%-5d\n" % (disp[0] + 1))
file.write("%20.16f %20.16f %20.16f\n" % (tuple(disp[1])))
# Subtract residual forces # Subtract residual forces
if not zero_forces==None: if zero_forces is not None:
forces[disp[2]] -= zero_forces forces[count] -= zero_forces
for f in forces[disp[2]]: for f in forces[count]:
file.write("%15.10f %15.10f %15.10f\n" % (tuple(f))) fp.write("%15.10f %15.10f %15.10f\n" % (tuple(f)))
# Show progress
if verbose:
print >> sys.stderr, "\n"
def mycmp(x, y): def mycmp(x, y):
return cmp(x[0], y[0]) return cmp(x[0], y[0])
# Sort by the atom numbering def parse_disp_yaml(filename="disp.yaml", return_cell=False):
# To remember the original order, the original index is added at 5th element.
def sort_displacements(displacements):
for i, disp in enumerate(displacements):
disp.append(i)
displacements.sort(mycmp)
return displacements
def parse_disp_yaml_with_supercell(filename='disp.yaml'):
try: try:
import yaml import yaml
except ImportError: except ImportError:
@ -269,45 +244,38 @@ def parse_disp_yaml_with_supercell(filename='disp.yaml'):
try: try:
from yaml import CLoader as Loader from yaml import CLoader as Loader
from yaml import CDumper as Dumper
except ImportError: except ImportError:
from yaml import Loader, Dumper from yaml import Loader
data = yaml.load(open(filename).read(), Loader=Loader) dataset = yaml.load(open(filename), Loader=Loader)
lattice = data['lattice'] natom = dataset['natom']
displacements = [] new_dataset = {}
for x in data['displacements']: new_dataset['natom'] = natom
displacements.append([x['atom']-1, x['displacement']]) new_first_atoms = []
positions = [x['position'] for x in data['atoms']] for first_atoms in dataset['displacements']:
symbols = [x['symbol'] for x in data['atoms']] first_atoms['atom'] -= 1
cell = Atoms(cell=lattice, atom1 = first_atoms['atom']
scaled_positions=positions, disp1 = first_atoms['displacement']
symbols=symbols, if 'direction' in first_atoms:
pbc=True) direct1 = first_atoms['direction']
new_first_atoms.append({'number': atom1,
'displacement': disp1,
'direction':direct1})
else:
new_first_atoms.append({'number': atom1, 'displacement': disp1})
new_dataset['first_atoms'] = new_first_atoms
return displacements, cell if return_cell:
lattice = dataset['lattice']
def parse_disp_yaml(filename='disp.yaml'): positions = [x['position'] for x in dataset['atoms']]
try: symbols = [x['symbol'] for x in dataset['atoms']]
import yaml cell = Atoms(cell=lattice,
except ImportError: scaled_positions=positions,
print "You need to install python-yaml." symbols=symbols,
exit(1) pbc=True)
return new_dataset, cell
try: else:
from yaml import CLoader as Loader return new_dataset
from yaml import CDumper as Dumper
except ImportError:
from yaml import Loader, Dumper
data = yaml.load(open(filename).read(), Loader=Loader)
natom = data['natom']
lattice = data['lattice']
displacements = []
for x in data['displacements']:
displacements.append([x['atom']-1, x['displacement']])
return displacements, natom
def parse_DISP(filename='DISP'): def parse_DISP(filename='DISP'):
disp = open(filename) disp = open(filename)

15
scripts/phonopy
View File

@ -411,29 +411,28 @@ if options.force_sets_mode or options.force_sets_zero_mode:
if file_exists('disp.yaml', log_level): if file_exists('disp.yaml', log_level):
if interface_mode == 'vasp': if interface_mode == 'vasp':
displacements, natom = file_IO.parse_disp_yaml('disp.yaml') displacements = file_IO.parse_disp_yaml('disp.yaml')
if interface_mode == 'wien2k': if interface_mode == 'wien2k':
displacements, supercell = file_IO.parse_disp_yaml_with_supercell('disp.yaml') displacements, supercell = file_IO.parse_disp_yaml(
'disp.yaml', return_cell=True)
for filename in args: for filename in args:
file_exists(filename, log_level) file_exists(filename, log_level)
num_disp_files = len(args) num_disp_files = len(args)
if options.force_sets_zero_mode: if options.force_sets_zero_mode:
num_disp_files -= 1 num_disp_files -= 1
if not len(displacements)== num_disp_files: if len(displacements['first_atoms']) != num_disp_files:
print_error_message("Number of files to be read don't match " print_error_message("Number of files to be read don't match "
"to number of displacements in disp.yaml.") "to number of displacements in disp.yaml.")
if log_level > 0: if log_level > 0:
print_end() print_end()
sys.exit(1) sys.exit(1)
if interface_mode == 'vasp': if interface_mode == 'vasp':
is_created = file_IO.write_FORCE_SETS_vasp( is_created = file_IO.write_FORCE_SETS_vasp(
args, args,
displacements, displacements,
natom,
'FORCE_SETS', 'FORCE_SETS',
is_zero_point=options.force_sets_zero_mode) is_zero_point=options.force_sets_zero_mode)
@ -450,7 +449,6 @@ if options.force_sets_mode or options.force_sets_zero_mode:
is_zero_point=options.force_sets_zero_mode, is_zero_point=options.force_sets_zero_mode,
is_distribute=(not options.is_wien2k_p1), is_distribute=(not options.is_wien2k_p1),
symprec=options.symprec) symprec=options.symprec)
if log_level > 0: if log_level > 0:
if is_created: if is_created:
@ -460,6 +458,9 @@ if options.force_sets_mode or options.force_sets_zero_mode:
print_end() print_end()
sys.exit(0) sys.exit(0)
##########################
# Read crystal structure #
##########################
# Parse VASP-type crystal structure # Parse VASP-type crystal structure
if interface_mode=='vasp': if interface_mode=='vasp':
if options.cell_poscar == None: if options.cell_poscar == None: