1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
"""
This module contains functionality for reading and writing an ASE
Atoms object in V_Sim 3.5+ ascii format.
"""
import numpy as np
from ase.utils import reader, writer
@reader
def read_v_sim(fd):
"""Import V_Sim input file.
Reads cell, atom positions, etc. from v_sim ascii file
"""
from ase import Atoms, units
from ase.geometry import cellpar_to_cell
import re
# Read comment:
fd.readline()
line = fd.readline() + ' ' + fd.readline()
box = line.split()
for i in range(len(box)):
box[i] = float(box[i])
keywords = []
positions = []
symbols = []
unit = 1.0
re_comment = re.compile(r'^\s*[#!]')
re_node = re.compile(r'^\s*\S+\s+\S+\s+\S+\s+\S+')
while True:
line = fd.readline()
if line == '':
break # EOF
p = re_comment.match(line)
if p is not None:
# remove comment character at the beginning of line
line = line[p.end():].replace(',', ' ').lower()
if line[:8] == "keyword:":
keywords.extend(line[8:].split())
elif re_node.match(line):
unit = 1.0
if not ("reduced" in keywords):
if (("bohr" in keywords) or ("bohrd0" in keywords) or
("atomic" in keywords) or ("atomicd0" in keywords)):
unit = units.Bohr
fields = line.split()
positions.append([unit * float(fields[0]),
unit * float(fields[1]),
unit * float(fields[2])])
symbols.append(fields[3])
if ("surface" in keywords) or ("freeBC" in keywords):
raise NotImplementedError
# create atoms object based on the information
if "angdeg" in keywords:
cell = cellpar_to_cell(box)
else:
unit = 1.0
if (("bohr" in keywords) or ("bohrd0" in keywords) or
("atomic" in keywords) or ("atomicd0" in keywords)):
unit = units.Bohr
cell = np.zeros((3, 3))
cell.flat[[0, 3, 4, 6, 7, 8]] = box[:6]
cell *= unit
if "reduced" in keywords:
atoms = Atoms(cell=cell, scaled_positions=positions)
else:
atoms = Atoms(cell=cell, positions=positions)
atoms.set_chemical_symbols(symbols)
return atoms
@writer
def write_v_sim(fd, atoms):
"""Write V_Sim input file.
Writes the atom positions and unit cell.
"""
from ase.geometry import cellpar_to_cell, cell_to_cellpar
# Convert the lattice vectors to triangular matrix by converting
# to and from a set of lengths and angles
cell = cellpar_to_cell(cell_to_cellpar(atoms.cell))
dxx = cell[0, 0]
dyx, dyy = cell[1, 0:2]
dzx, dzy, dzz = cell[2, 0:3]
fd.write('===== v_sim input file created using the'
' Atomic Simulation Environment (ASE) ====\n')
fd.write('{0} {1} {2}\n'.format(dxx, dyx, dyy))
fd.write('{0} {1} {2}\n'.format(dzx, dzy, dzz))
# Use v_sim 3.5 keywords to indicate scaled positions, etc.
fd.write('#keyword: reduced\n')
fd.write('#keyword: angstroem\n')
if np.alltrue(atoms.pbc):
fd.write('#keyword: periodic\n')
elif not np.any(atoms.pbc):
fd.write('#keyword: freeBC\n')
elif np.array_equiv(atoms.pbc, [True, False, True]):
fd.write('#keyword: surface\n')
else:
raise Exception(
'Only supported boundary conditions are full PBC,'
' no periodic boundary, and surface which is free in y direction'
' (i.e. Atoms.pbc = [True, False, True]).')
# Add atoms (scaled positions)
for position, symbol in zip(atoms.get_scaled_positions(),
atoms.get_chemical_symbols()):
fd.write('{0} {1} {2} {3}\n'.format(
position[0], position[1], position[2], symbol))
|
