#
# The high-throughput toolkit (httk)
# Copyright (C) 2012-2013 Rickard Armiento
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import time
try:
import cfractions as fractions
except ImportError:
import fractions
from httk.core import FracVector, MutableFracVector
from httk.core.basic import is_sequence
from httk.atomistic import spacegrouputils
[docs]
def sort_coordgroups(coordgroups, individual_data):
counts = [len(x) for x in coordgroups]
newcoordgroups = []
newindividual_data = []
for group in range(len(counts)):
order = sorted(range(counts[group]), key=lambda x: (coordgroups[group][x][0], coordgroups[group][x][1], coordgroups[group][x][2]))
newcoordgroups.append(coordgroups[group][order])
if individual_data is not None:
newindividual_data.append([individual_data[group][i] for i in order])
if individual_data is None:
return coordgroups.stack(newcoordgroups), newindividual_data
else:
return coordgroups.stack(newcoordgroups), None
[docs]
def coords_and_occupancies_to_coordgroups_and_assignments(coords, occupancies):
if len(occupancies) != len(coords):
raise Exception("Occupations and coords needs to be of the same length.")
if len(occupancies) == 0:
return [], FracVector((), 1)
coordgroups = []
group_occupancies = []
for i in range(len(occupancies)):
try:
idx = group_occupancies.index(occupancies[i])
except ValueError:
idx = len(group_occupancies)
group_occupancies.append(occupancies[i])
coordgroups.append([])
coordgroups[idx].append(coords[i])
return coordgroups, group_occupancies
[docs]
def coords_to_coordgroups(coords, counts):
coordgroups = []
idx = 0
for count in counts:
coordgroups.append(coords[idx:count+idx])
idx += count
return coordgroups
[docs]
def coordgroups_to_coords(coordgroups):
coords = []
counts = [len(x) for x in coordgroups]
for group in coordgroups:
for i in range(len(group)):
coords.append(group[i])
coords = FracVector.stack_vecs(coords)
return coords, counts
[docs]
def coordgroups_cartesian_to_reduced(coordgroups, basis):
basis = FracVector.use(basis)
cellinv = basis.inv()
newcoordgroups = []
for coordgroup in coordgroups:
newcoordgroup = coordgroup*cellinv
newcoordgroups.append(newcoordgroup)
return newcoordgroups
[docs]
def coords_and_counts_to_coordgroups(coords, counts):
coordgroups = []
idx = 0
for count in counts:
coordgroups.append(coords[idx:count+idx])
idx += count
return coordgroups
[docs]
def coordswap(fromidx, toidx, cell, coordgroups):
new_coordgroups = []
for group in coordgroups:
coords = MutableFracVector.from_FracVector(group)
rows = coords[:, toidx]
coords[:, toidx] = coords[:, fromidx]
coords[:, fromidx] = rows
new_coordgroups.append(coords.to_FracVector())
coordgroups = FracVector.create(new_coordgroups)
cell = MutableFracVector.from_FracVector(cell)
row = cell[toidx]
cell[toidx] = cell[fromidx]
cell[fromidx] = row
cell = cell.to_FracVector()
return (cell, coordgroups)
[docs]
def clean_coordgroups_and_assignments(coordgroups, assignments):
if len(assignments) != len(coordgroups):
raise Exception("Occupations and coords needs to be of the same length.")
if len(assignments) == 0:
return [], FracVector((), 1)
new_coordgroups = []
new_assignments = []
for i in range(len(assignments)):
for j in range(len(new_assignments)):
if assignments[i] == new_assignments[j]:
idx = j
new_coordgroups[idx] = FracVector.chain_vecs([new_coordgroups[idx], coordgroups[i]])
break
else:
new_coordgroups.append(coordgroups[i])
new_assignments.append(assignments[i])
idx = len(coordgroups)-1
return new_coordgroups, new_assignments
# TODO: Cleanup formula generation
[docs]
def abstract_symbol(count):
if count < 27:
return chr(64+count)
my = (count-1) % 26
rec = (count-1) // 26
return abstract_symbol(rec)+chr(97+my)
[docs]
def coords_reduced_to_cartesian(cell, coords):
cell = FracVector.use(cell)
newcoords = coords*cell
return newcoords
[docs]
def coordgroups_reduced_to_cartesian(cell, coordgroups):
cell = FracVector.use(cell)
newcoordgroups = []
for coordgroup in coordgroups:
newcoordgroup = coordgroup*cell
newcoordgroups.append(newcoordgroup)
return newcoordgroups
[docs]
def periodicity_to_pbc(periodicity):
if is_sequence(periodicity):
pbc = [bool(x) for x in periodicity]
else:
pbc = [False]*periodicity + [True]*(3-periodicity)
return pbc
[docs]
def pbc_to_nonperiodic_vecs(pbc):
nonperiodic_vecs = sum(pbc)
if (pbc[0] and not (pbc[1] or pbc[2])) or (pbc[1] and not pbc[2]):
raise Exception("cellutils pbc_to_nonperiodic_vecs: cannot represent as nonperiodic vecs, need non-periodic vecs before periodic ones.")
return nonperiodic_vecs
[docs]
def structure_reduced_coordgroups_to_representative(coordgroups, cell,
spacegroup, backends=['isotropy', 'spglib']):
for backend in backends:
if backend == 'isotropy':
try:
from httk.external import isotropy_ext
return isotropy_ext.uc_reduced_coordgroups_process_with_isotropy(coordgroups, cell, spacegroup, get_wyckoff=True)
except ImportError:
raise
pass
if backend == 'spglib':
try:
from httk.external import pyspglib_ext
return pyspglib_ext.uc_reduced_coordgroups_process_with_isotropy(coordgroups, cell, spacegroup, get_wyckoff=True)
except ImportError:
raise
pass
raise Exception("structure_reduced_coordgroups_to_representative: None of the available backends available.")
[docs]
def sites_tidy(sites, backends=['platon']):
for backend in backends:
if backend == 'platon':
try:
from httk.external import platon_ext
return platon_ext.sites_tidy(sites)
except ImportError:
raise
pass
raise Exception("structure_tidy: None of the available backends available.")
[docs]
def coordgroups_reduced_to_unitcell(coordgroups, hall_symbol, eps=fractions.Fraction(1, 1000)):
symops = spacegrouputils.get_symops(hall_symbol)
# Number of grid subdivisions along each axis, e.g. if eps=1/1000 => g=1000
g = max(10, int(1 / float(eps)))
def neighbor_cells(i, j, k):
for di in (-1, 0, 1):
for dj in (-1, 0, 1):
for dk in (-1, 0, 1):
yield ((i + di) % g, (j + dj) % g, (k + dk) % g)
newcoordgroups = []
for coordgroup in coordgroups:
newcoordgroup = []
# We'll keep one grid-hash dictionary per coordgroup
# Key: (i,j,k) cell index -> Value: list of accepted coords in that cell
cell_dict = {}
for symop in symops:
rotcoords = coordgroup * symop[0].T()
for coord in rotcoords:
finalcoord = (coord + symop[1]).normalize()
i, j, k = int(finalcoord[0] * g) % g, int(finalcoord[1] * g) % g, int(finalcoord[2] * g) % g
too_close = False
for (ni, nj, nk) in neighbor_cells(i, j, k):
for c in cell_dict.get((ni, nj, nk), []):
if (c - finalcoord).normalize_half().lengthsqr() < eps:
too_close = True
break
if too_close:
break
if not too_close:
newcoordgroup.append(finalcoord)
cell_dict.setdefault((i, j, k), []).append(finalcoord)
newcoordgroups.append(newcoordgroup)
return newcoordgroups
[docs]
def coordgroups_reduced_to_unitcell_old(coordgroups, hall_symbol, eps=fractions.Fraction(1,1000)):
# TODO: This routine is now a bit faster, but maybe can be accelerated further
# One may also be concerned about how much memory it uses now.
symops = spacegrouputils.get_symops(hall_symbol)
#print("ENTER CRDTU",sum(len(x) for x in coordgroups),len(symops))
#start = time.process_time()
#count = 0
newcoordgroups = []
for coordgroup in coordgroups:
newcoordgroup = set()
for symop in symops:
rotcoords = coordgroup*(symop[0].T())
for coord in rotcoords:
#count += 1
finalcoord = (coord+symop[1]).normalize()
newcoordgroup.add(finalcoord)
if finalcoord not in newcoordgroup:
for checkcoord in newcoordgroup:
if (checkcoord-finalcoord).normalize_half().lengthsqr() < eps:
break
else:
newcoordgroup.add(finalcoord)
#newcoordgroup = sorted(list(newcoordgroup), key=lambda x: (x[0], x[1], x[2]))
#newcoordgroup2 = [newcoordgroup[0]]
#lastcoord = newcoordgroup[0]
#for coord in newcoordgroup[1:]:
# if (coord-lastcoord).normalize_half().lengthsqr() >= eps:
# newcoordgroup2.append(coord)
# lastcoord = coord
newcoordgroups += [newcoordgroup]
#old = coordgroups_reduced_to_unitcell_old(coordgroups, hall_symbol, eps)
#if newcoordgroups != old:
# print("\n\nONE\n",FracVector.create(old).to_floats())
# print("\n\nTWO\n",FracVector.create(newcoordgroups).to_floats())
# raise Exception("DIDN'T WORK")
return newcoordgroups
[docs]
def main():
cell = FracVector.create([[1, 1, 0], [1, 0, 1], [0, 1, 1]])
coordgroups = FracVector.create([[[2, 3, 5], [3, 5, 4]], [[4, 6, 7]]])
assignments = [2, 5]
print(cell, coordgroups)
cell, coordgroups = coordswap(0, 2, cell, coordgroups)
print(cell, coordgroups)
pass
if __name__ == "__main__":
main()