Commit 167bb7c6 authored by Carl Poelking's avatar Carl Poelking
Browse files

Connectivity tables etc.

parent e8b31758
......@@ -166,9 +166,11 @@ class KernelAdaptorSpecificUnique(object):
IX = np.zeros((0,0), dtype='float64') # feature matrix
dimX = -1
IR = np.zeros((0,0), dtype='float64') # position matrix
types = []
for atomic_i in spectrum:
Xi_unnorm, Xi_norm = self.adaptScalar(atomic_i)
Ri = atomic_i.getCenter().pos
types.append(atomic_i.getCenter().type)
dimX = Xi_norm.shape[0]
if not IX.any():
IX = np.copy(Xi_norm) # TODO Is this necessary?
......@@ -182,7 +184,7 @@ class KernelAdaptorSpecificUnique(object):
IR.resize((i+1, 3))
IR[-1,:] = Ri
if return_pos_matrix:
return IX, IR
return IX, IR, types
else:
return IX
def adaptScalar(self, atomic):
......
......@@ -42,6 +42,24 @@ def mp_compute_column_block(gi, gj_list, kfct):
krow.append(k)
return krow
def compute_upper_triangle(
kfct,
g_list,
**kwargs):
dim = len(g_list)
kmat = np.zeros((dim,dim), dtype='float64')
kfct_primed = fct.partial(
kfct,
**kwargs)
for i in range(dim):
gi = g_list[i]
for j in range(i, dim):
gj = g_list[j]
kij = kfct(gi, gj)
kmat[i,j] = kij
kmat[j,i] = kij
return kmat
def mp_compute_upper_triangle(
kfct,
g_list,
......
......@@ -20,76 +20,127 @@ def structures_from_xyz(xyz_file, do_partition=True, add_fragment_com=True):
do_partition=do_partition,
add_fragment_com=add_fragment_com)
def structures_from_ase(ase_configs, do_partition=True, add_fragment_com=True, use_center_of_geom=False, log=None):
# TODO def structures_from_ase in addition to structure_from_ase
def structure_from_ase(
config,
do_partition=True,
add_fragment_com=True,
use_center_of_geom=False,
log=None):
# NOTE Center of mass is computed without considering PBC => Requires unwrapped coordinates
# Create structures
structures = []
for config in ase_configs:
# System properties
label = config.info['label']
positions = config.get_positions()
types = config.get_chemical_symbols()
if log: log << log.back << "Reading '%s'" % label << log.flush
# Simulation cell
if config.pbc.all():
box = np.array([ase_config.cell[0], ase_config.cell[1], ase_config.cell[2]])
elif not config.pbc.any():
box = np.zeros((3,3))
else:
raise NotImplementedError("<structures_from_xyz> Partial periodicity not implemented.")
# Partition
if do_partition:
frags, top = partition.PartitionStructure(types, positions)
positions = [ atm.xyz for atm in itertools.chain(*frags) ]
types = [ atm.e for atm in itertools.chain(*frags) ]
else:
top = [('SEG', 1, len(positions))]
# Check particle count consistent with top
atom_count = 0
for section in top:
atom_count += section[1]*section[2]
assert atom_count == len(positions) # Does topology match structure?
# Create segments, particles
structure = soap.Structure(label)
structure.box = box
atom_idx = 0
for section in top:
seg_type = section[0]
n_segs = section[1]
n_atoms = section[2]
for i in range(n_segs):
structure = None
frag_bond_matrix = None
atom_bond_matrix = None
frag_labels = []
atom_labels = []
# System properties
label = config.info['label']
positions = config.get_positions()
types = config.get_chemical_symbols()
if log: log << log.back << "Reading '%s'" % label << log.flush
# Simulation cell
if config.pbc.all():
box = np.array([config.cell[0], config.cell[1], config.cell[2]])
elif not config.pbc.any():
box = np.zeros((3,3))
else:
raise NotImplementedError("<structures_from_xyz> Partial periodicity not implemented.")
# Partition
if do_partition:
# Partition => Frags, top, reordered positions
frags, top = partition.PartitionStructure(types, positions, outfile_gro='%s.gro' % label)
atms = [ atm for atm in itertools.chain(*frags) ]
positions = [ atm.xyz for atm in itertools.chain(*frags) ]
types = [ atm.e for atm in itertools.chain(*frags) ]
# Generate fragment labels
for frag in frags:
label = '-'.join(atm.e for atm in frag)
frag_labels.append(label)
# Reorder map
id_reorder_map = {}
for atm in atms:
id_reorder_map[atm.id_initial] = atm.id
# Connectivity matrix: fragments
frag_bond_matrix = np.zeros((len(frags),len(frags)), dtype=bool)
for i in range(len(frags)):
frag_bond_matrix[i,i] = True
for j in range(i+1, len(frags)):
frags_are_bonded = False
for ai in frags[i]:
for aj in frags[j]:
if aj in ai.bonded:
frags_are_bonded = True
break
if frags_are_bonded: break
frag_bond_matrix[i,j] = frags_are_bonded
frag_bond_matrix[j,i] = frags_are_bonded
# Connectivity matrix: atoms
atom_bond_matrix = np.zeros((len(positions),len(positions)), dtype=bool)
for i in range(len(atms)):
atom_bond_matrix[i,i] = True
ai = atms[i]
for j in range(i+1, len(atms)):
atoms_are_bonded = False
aj = atms[j]
if aj in ai.bonded:
atoms_are_bonded = True
atom_bond_matrix[i,j] = atoms_are_bonded
atom_bond_matrix[j,i] = atoms_are_bonded
# Reorder ASE atoms
config = ase.Atoms(
sorted(config, key = lambda atm: id_reorder_map[atm.index+1]),
info=config.info,
cell=config.cell,
pbc=config.pbc)
else:
top = [('SEG', 1, len(positions))]
# Check particle count consistent with top
atom_count = 0
for section in top:
atom_count += section[1]*section[2]
assert atom_count == len(positions) # Does topology match structure?
# Create segments, particles
structure = soap.Structure(label)
structure.box = box
atom_idx = 0
for section in top:
seg_type = section[0]
n_segs = section[1]
n_atoms = section[2]
for i in range(n_segs):
segment = structure.addSegment()
segment.name = seg_type
# Add particles, compute CoM
com = np.array([0.,0.,0.])
com_weight_total = 0
for j in range(n_atoms):
particle = structure.addParticle(segment)
particle.pos = positions[atom_idx]
particle.weight = 1.
particle.sigma = 0.5
particle.type = types[atom_idx]
atom_labels.append(particle.type)
# Compute CoMass/CoGeom
if use_center_of_geom:
com_weight = 1.
else:
com_weight = elements.periodic_table[particle.type].mass
com_weight_total += com_weight
com = com + com_weight*positions[atom_idx]
atom_idx += 1
com = com/com_weight_total
# Add CoM particle if requested
if add_fragment_com:
segment = structure.addSegment()
segment.name = seg_type
# Add particles, compute CoM
com = np.array([0.,0.,0.])
com_weight_total = 0
for j in range(n_atoms):
particle = structure.addParticle(segment)
particle.pos = positions[atom_idx]
particle.weight = 1.
particle.sigma = 0.5
particle.type = types[atom_idx]
# Compute CoMass/CoGeom
if use_center_of_geom:
com_weight = 1.
else:
com_weight = elements.periodic_table[particle.type].mass
com_weight_total += com_weight
com = com + com_weight*positions[atom_idx]
atom_idx += 1
com = com/com_weight_total
# Add CoM particle if requested
if add_fragment_com:
segment = structure.addSegment()
segment.name = "%s.COM" % seg_type
particle = structure.addParticle(segment)
particle.pos = com
particle.weight = 0.
particle.sigma = 0.5
particle.type = "COM"
structures.append(structure)
segment.name = "%s.COM" % seg_type
particle = structure.addParticle(segment)
particle.pos = com
particle.weight = 0.
particle.sigma = 0.5
particle.type = "COM"
if log: log << log.endl
return structures
return config, structure, top, frag_bond_matrix, atom_bond_matrix, frag_labels, atom_labels
def write_xyz(xyz_file, structure):
ofs = open(xyz_file, 'w')
......
......@@ -79,6 +79,7 @@ class AtomXyz(object):
def __init__(self, e, xyz, id):
# PROPERTIES
self.e = e
self.id_initial = id
self.id = id
self.xyz = xyz
self.x = xyz[0]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment