#!/usr/bin/env python
# PE, 2015/06/04
description_string =
"""The script reads an ASE compatible input file containing one or
more configurations of positions and forces, e.g. a VASP OUTCAR file
or a .traj file. It then writes one or more corresponding files in a
POSCAR like format including forces, which is suitable for force
matching using the atomicrex code.
"""
from ase import *
from ase.io import *
#------------------------------------------------
# Handle the command line
import argparse
parser = argparse.ArgumentParser(description=description_string)
parser.add_argument('filename', help='input structure in ASE readable format; the file can contain multiple configurations')
parser.add_argument('-o', '--outfile', default='FPOSCAR.new', help='name of output file')
args = parser.parse_args()
#------------------------------------------------
def write_configuration(outfile, conf, direct=False, comment='', form='%21.16f'):
with open(outfile, 'w') as f:
# Element names
elements = []
for elem in conf.get_chemical_symbols():
if not elem in elements: elements.append(elem)
# Comment line (here we store the energy per atom)
energy = conf.get_potential_energy()
energy /= conf.get_number_of_atoms()
f.write('%.6g\n' % energy)
# Cell metric
f.write('1.0\n')
for vec in conf.cell:
for v in vec: f.write((' '+form) % v)
f.write('\n')
# Element names
for elem in elements: f.write(' %2s' % elem)
f.write('\n')
# Number of elements per type
for elem in elements:
f.write(' %d' % conf.get_chemical_symbols().count(elem))
f.write('\n')
# Scaled or cartesian coordinates
if direct:
coordinates = atoms.get_scaled_positions()
f.write('Direct\n')
else:
coordinates = conf.get_positions()
f.write('Cartesian\n')
# Coordinates and forces
for coord,force in zip(coordinates, conf.get_forces()):
for r in coord: f.write((' '+form) % r)
f.write(' ')
for r in force: f.write((' '+form) % r)
f.write('\n')
#------------------------------------------------
# Read and write configuration
fname = args.filename
try:
confs = read(fname, ':')
except:
print "ERROR: file %s could not be read"%fname
raise
for k,conf in enumerate(confs):
if len(confs) > 1:
outfile = args.outfile+'_%d' % k
else:
outfile = args.outfile
write_configuration(outfile, conf)