Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def read_geometry(filename, dir='.'):
fh = open(os.path.join(dir, filename), 'rb')
lines = list(filter(read_geometry_filter, fh.readlines()))
# return geometry in ASE format
geometry = []
for line in lines:
sline = line.split()
# find chemical symbol (the symbols in the file are lowercase)
symbol = sline[-1]
for s in chemical_symbols:
if symbol == s.lower():
symbol = s
break
geometry.append(Atom(symbol=symbol, position=sline[:-1]))
fh.close()
atoms = Atoms(geometry)
atoms.set_positions(atoms.get_positions()*Bohr) # convert to Angstrom
return atoms
table[n] = ':' + table[n] + ':'
table = '\n'.join([s for s in table])
# split into compounds
# http://simonwillison.net/2003/Oct/26/reSplit/
# http://stackoverflow.com/questions/647655/python-regex-split-and-special-character
table = re.compile('(:.*:)').split(table)
# remove empty elements
table = [l.strip() for l in table]
table = [l for l in table if len(l) > 1]
# extract compounds
for n in range(0, len(table), 2):
compound = table[n].replace(':', '').replace(' ', '_')
geometry = []
for atom in table[n+1].split('\n'):
geometry.append(Atom(symbol=atom.split()[0], position=atom.split()[1:]))
atoms = Atoms(geometry)
# set the charge and magnetic moment on the heaviest atom (better ideas?)
heaviest = max([a.get_atomic_number() for a in atoms])
heaviest_index = [a.get_atomic_number() for a in atoms].index(heaviest)
charge = 0.0
if abs(charge) > 0.0:
charges = [0.0 for a in atoms]
charges[heaviest_index] = charge
atoms.set_initial_charges(charges)
if compound in [ # see corresponding articles
'Ti(BH4)3', # TM1R2006
'V(NMe2)4', # TM1R2006
'Cu(acac)2', # TM1R2006
'Nb(Cp)(C7H7)_Cs', # TM2R2007
'CdMe_C3v', # TM2R2007
]:
multiplicity = 2.0
# Now we have to account for the possibility of there being CASTEP
# 'custom' species amongst the symbols
custom_species = None
for a in data_dict['atoms']['atom']:
spec_custom = a['species'].split(':', 1)
if len(spec_custom) > 1 and custom_species is None:
# Add it to the custom info!
custom_species = list(symbols)
symbols.append(spec_custom[0])
positions.append(a['position'])
indices.append(a['index'])
labels.append(a['label'])
if custom_species is not None:
custom_species.append(a['species'])
atoms = Atoms(cell=cell,
pbc=pbc,
symbols=symbols,
positions=positions)
# Add custom species if present
if custom_species is not None:
atoms.new_array('castep_custom_species', np.array(custom_species))
# Add the spacegroup, if present and recognizable
if 'symmetry' in data_dict['atoms']:
try:
spg = Spacegroup(data_dict['atoms']['symmetry'][0])
except SpacegroupNotFoundError:
# Not found
spg = Spacegroup(1) # Most generic one
atoms.info['spacegroup'] = spg
def unpack_reftraj_str_to_atoms(data):
lines = data.split(b'\n')
label = int(lines[0])
n_atoms = int(lines[1])
at = Atoms(symbols=[' ']*n_atoms, cell=np.eye(3))
at.info['label'] = label
for i in range(3):
at.cell[:, i] = [float(x) for x in lines[i].split()]
for i, line in enumerate(lines[4:]):
t = [float(x) for x in line.split()]
at.positions[i, :] = np.dot(t, at.cell)
return at
def build_molecule(self, name):
args = self.args
try:
# Known molecule or atom?
atoms = molecule(name)
except NotImplementedError:
symbols = string2symbols(name)
if len(symbols) == 1:
Z = atomic_numbers[symbols[0]]
magmom = ground_state_magnetic_moments[Z]
atoms = Atoms(name, magmoms=[magmom])
elif len(symbols) == 2:
# Dimer
if args.bond_length is None:
b = (covalent_radii[atomic_numbers[symbols[0]]] +
covalent_radii[atomic_numbers[symbols[1]]])
else:
b = args.bond_length
atoms = Atoms(name, positions=[(0, 0, 0),
(b, 0, 0)])
else:
raise ValueError('Unknown molecule: ' + name)
else:
if len(atoms) == 2 and args.bond_length is not None:
atoms.set_distance(0, 1, args.bond_length)
if args.unit_cell is None:
def colored(self, elements):
res = Atoms()
res.set_cell(self.get_cell())
for atom in self:
elem = self.types[atom.tag]
if elem in elements:
elem = elements[elem]
res.append(Atom(elem, atom.position))
return res
"""
b = sqrt(3) * C_C / 4
arm_unit = Atoms(main_element + '4',
pbc=(1, 0, 1),
cell=[4 * b, 0, 3 * C_C])
arm_unit.positions = [[0, 0, 0],
[b * 2, 0, C_C / 2.],
[b * 2, 0, 3 * C_C / 2.],
[0, 0, 2 * C_C]]
zz_unit = Atoms(main_element + '2',
pbc=(1, 0, 1),
cell=[3 * C_C / 2.0, 0, b * 4])
zz_unit.positions = [[0, 0, 0],
[C_C / 2.0, 0, b * 2]]
atoms = Atoms()
if type == 'zigzag':
edge_index0 = np.arange(m) * 2
edge_index1 = (n - 1) * m * 2 + np.arange(m) * 2 + 1
if magnetic:
mms = np.zeros(m * n * 2)
for i in edge_index0:
mms[i] = initial_mag
for i in edge_index1:
mms[i] = -initial_mag
for i in range(n):
layer = zz_unit.repeat((1, 1, m))
layer.positions[:, 0] += 3 * C_C / 2 * i
if i % 2 == 1:
def read_gromacs(filename):
""" From:
http://manual.gromacs.org/current/online/gro.html
C format
"%5d%-5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f"
python: starting from 0, including first excluding last
0:4 5:10 10:15 15:20 20:28 28:36 36:44 44:52 52:60 60:68
Import gromacs geometry type files (.gro).
Reads atom positions,
velocities(if present) and
simulation cell (if present)
"""
atoms = Atoms()
filed = open(filename, 'r')
lines = filed.readlines()
filed.close()
positions = []
gromacs_velocities = []
symbols = []
tags = []
gromacs_residuenumbers = []
gromacs_residuenames = []
gromacs_atomtypes = []
sym2tag = {}
tag = 0
for line in (lines[2:-1]):
#print line[0:5]+':'+line[5:11]+':'+line[11:15]+':'+line[15:20]
# it is not a good idea to use the split method with gromacs input
# since the fields are defined by a fixed column number. Therefore,
def _write(filename, fd, format, io, images, parallel=None, append=False,
**kwargs):
if isinstance(images, Atoms):
images = [images]
if io.single:
if len(images) > 1:
raise ValueError('{}-format can only store 1 Atoms object.'
.format(format))
images = images[0]
if io.write is None:
raise ValueError("Can't write to {}-format".format(format))
# Special case for json-format:
if format == 'json' and (len(images) > 1 or append):
if filename is not None:
io.write(filename, images, append=append, **kwargs)
return
[0.188154, 2.210362, 0.883614],
[0.188154, 2.210362, -0.883614],
[-0.188154, -2.210362, 0.883614],
[-0.188154, -2.210362, -0.883614],
[1.247707, -0.07266, -0.877569],
[1.247707, -0.07266, 0.877569],
[-1.247707, 0.07266, -0.877569],
[-1.247707, 0.07266, 0.877569]],
'symbols': 'CCCCHHHHHHHHHH',
'thermal correction': 4.2633}}
# all constituent atoms
atoms_g22 = []
for f in data.keys():
s = Atoms(symbols=data[f]['symbols'], positions=data[f]['positions'])
for a in s:
atoms_g22.append(a.symbol)
# unique atoms
atoms_g22 = list(set(atoms_g22))
# add remaining atoms from G2_1
from ase.data.g2_1 import data as data1
for a in atoms_g22:
if not a in data.keys():
data[a] = data1[a]