Coverage for /builds/kinetik161/ase/ase/io/xsd.py: 69.28%
166 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import xml.etree.ElementTree as ET
2from xml.dom import minidom
4import numpy as np
6from ase import Atoms
7from ase.utils import writer
10def read_xsd(fd):
11 tree = ET.parse(fd)
12 root = tree.getroot()
14 atomtreeroot = root.find('AtomisticTreeRoot')
15 # if periodic system
16 if atomtreeroot.find('SymmetrySystem') is not None:
17 symmetrysystem = atomtreeroot.find('SymmetrySystem')
18 mappingset = symmetrysystem.find('MappingSet')
19 mappingfamily = mappingset.find('MappingFamily')
20 system = mappingfamily.find('IdentityMapping')
22 coords = []
23 cell = []
24 formula = ''
26 for atom in system:
27 if atom.tag == 'Atom3d':
28 symbol = atom.get('Components')
29 formula += symbol
31 xyz = atom.get('XYZ')
32 if xyz:
33 coord = [float(coord) for coord in xyz.split(',')]
34 else:
35 coord = [0.0, 0.0, 0.0]
36 coords.append(coord)
37 elif atom.tag == 'SpaceGroup':
38 avec = [float(vec) for vec in atom.get('AVector').split(',')]
39 bvec = [float(vec) for vec in atom.get('BVector').split(',')]
40 cvec = [float(vec) for vec in atom.get('CVector').split(',')]
42 cell.append(avec)
43 cell.append(bvec)
44 cell.append(cvec)
46 atoms = Atoms(formula, cell=cell, pbc=True)
47 atoms.set_scaled_positions(coords)
48 return atoms
49 # if non-periodic system
50 elif atomtreeroot.find('Molecule') is not None:
51 system = atomtreeroot.find('Molecule')
53 coords = []
54 formula = ''
56 for atom in system:
57 if atom.tag == 'Atom3d':
58 symbol = atom.get('Components')
59 formula += symbol
61 xyz = atom.get('XYZ')
62 coord = [float(coord) for coord in xyz.split(',')]
63 coords.append(coord)
65 atoms = Atoms(formula, pbc=False)
66 atoms.set_scaled_positions(coords)
67 return atoms
70def CPK_or_BnS(element):
71 """Determine how atom is visualized"""
72 if element in ['C', 'H', 'O', 'S', 'N']:
73 visualization_choice = 'Ball and Stick'
74 else:
75 visualization_choice = 'CPK'
76 return visualization_choice
79def SetChild(parent, childname, props):
80 Child = ET.SubElement(parent, childname)
81 for key in props:
82 Child.set(key, props[key])
83 return Child
86def SetBasicChilds():
87 """
88 Basic property setup for Material Studio File
89 """
90 XSD = ET.Element('XSD')
91 XSD.set('Version', '6.0')
93 ATR = SetChild(XSD, 'AtomisticTreeRoot', dict(
94 ID='1',
95 NumProperties='40',
96 NumChildren='1',
97 ))
98 SetChild(ATR, 'Property', dict(
99 DefinedOn='ClassicalEnergyHolder',
100 Name='AngleEnergy',
101 Type='Double',
102 ))
103 SetChild(ATR, 'Property', dict(
104 DefinedOn='ClassicalEnergyHolder',
105 Name='BendBendEnergy',
106 Type='Double',
107 ))
108 SetChild(ATR, 'Property', dict(
109 DefinedOn='ClassicalEnergyHolder',
110 Name='BendTorsionBendEnergy',
111 Type='Double',
112 ))
113 SetChild(ATR, 'Property', dict(
114 DefinedOn='ClassicalEnergyHolder',
115 Name='BondEnergy',
116 Type='Double',
117 ))
118 SetChild(ATR, 'Property', dict(
119 DefinedOn='Atom',
120 Name='EFGAsymmetry',
121 Type='Double',
122 ))
123 SetChild(ATR, 'Property', dict(
124 DefinedOn='Atom',
125 Name='EFGQuadrupolarCoupling',
126 Type='Double',
127 ))
128 SetChild(ATR, 'Property', dict(
129 DefinedOn='ClassicalEnergyHolder',
130 Name='ElectrostaticEnergy',
131 Type='Double',
132 ))
133 SetChild(ATR, 'Property', dict(
134 DefinedOn='GrowthFace',
135 Name='FaceMillerIndex',
136 Type='MillerIndex',
137 ))
138 SetChild(ATR, 'Property', dict(
139 DefinedOn='GrowthFace',
140 Name='FacetTransparency',
141 Type='Float',
142 ))
143 SetChild(ATR, 'Property', dict(
144 DefinedOn='Bondable',
145 Name='Force',
146 Type='CoDirection',
147 ))
148 SetChild(ATR, 'Property', dict(
149 DefinedOn='ClassicalEnergyHolder',
150 Name='HydrogenBondEnergy',
151 Type='Double',
152 ))
153 SetChild(ATR, 'Property', dict(
154 DefinedOn='Bondable',
155 Name='ImportOrder',
156 Type='UnsignedInteger',
157 ))
158 SetChild(ATR, 'Property', dict(
159 DefinedOn='ClassicalEnergyHolder',
160 Name='InversionEnergy',
161 Type='Double',
162 ))
163 SetChild(ATR, 'Property', dict(
164 DefinedOn='Atom',
165 Name='IsBackboneAtom',
166 Type='Boolean',
167 ))
168 SetChild(ATR, 'Property', dict(
169 DefinedOn='Atom',
170 Name='IsChiralCenter',
171 Type='Boolean',
172 ))
173 SetChild(ATR, 'Property', dict(
174 DefinedOn='Atom',
175 Name='IsOutOfPlane',
176 Type='Boolean',
177 ))
178 SetChild(ATR, 'Property', dict(
179 DefinedOn='BestFitLineMonitor',
180 Name='LineExtentPadding',
181 Type='Double',
182 ))
183 SetChild(ATR, 'Property', dict(
184 DefinedOn='Linkage',
185 Name='LinkageGroupName',
186 Type='String',
187 ))
188 SetChild(ATR, 'Property', dict(
189 DefinedOn='PropertyList',
190 Name='ListIdentifier',
191 Type='String',
192 ))
193 SetChild(ATR, 'Property', dict(
194 DefinedOn='Atom',
195 Name='NMRShielding',
196 Type='Double',
197 ))
198 SetChild(ATR, 'Property', dict(
199 DefinedOn='ClassicalEnergyHolder',
200 Name='NonBondEnergy',
201 Type='Double',
202 ))
203 SetChild(ATR, 'Property', dict(
204 DefinedOn='Bondable',
205 Name='NormalMode',
206 Type='Direction',
207 ))
208 SetChild(ATR, 'Property', dict(
209 DefinedOn='Bondable',
210 Name='NormalModeFrequency',
211 Type='Double',
212 ))
213 SetChild(ATR, 'Property', dict(
214 DefinedOn='Bondable',
215 Name='OrbitalCutoffRadius',
216 Type='Double',
217 ))
218 SetChild(ATR, 'Property', dict(
219 DefinedOn='BestFitPlaneMonitor',
220 Name='PlaneExtentPadding',
221 Type='Double',
222 ))
223 SetChild(ATR, 'Property', dict(
224 DefinedOn='ClassicalEnergyHolder',
225 Name='PotentialEnergy',
226 Type='Double',
227 ))
228 SetChild(ATR, 'Property', dict(
229 DefinedOn='ScalarFieldBase',
230 Name='QuantizationValue',
231 Type='Double',
232 ))
233 SetChild(ATR, 'Property', dict(
234 DefinedOn='ClassicalEnergyHolder',
235 Name='RestraintEnergy',
236 Type='Double',
237 ))
238 SetChild(ATR, 'Property', dict(
239 DefinedOn='ClassicalEnergyHolder',
240 Name='SeparatedStretchStretchEnergy',
241 Type='Double',
242 ))
243 SetChild(ATR, 'Property', dict(
244 DefinedOn='Trajectory',
245 Name='SimulationStep',
246 Type='Integer',
247 ))
248 SetChild(ATR, 'Property', dict(
249 DefinedOn='ClassicalEnergyHolder',
250 Name='StretchBendStretchEnergy',
251 Type='Double',
252 ))
253 SetChild(ATR, 'Property', dict(
254 DefinedOn='ClassicalEnergyHolder',
255 Name='StretchStretchEnergy',
256 Type='Double',
257 ))
258 SetChild(ATR, 'Property', dict(
259 DefinedOn='ClassicalEnergyHolder',
260 Name='StretchTorsionStretchEnergy',
261 Type='Double',
262 ))
263 SetChild(ATR, 'Property', dict(
264 DefinedOn='ClassicalEnergyHolder',
265 Name='TorsionBendBendEnergy',
266 Type='Double',
267 ))
268 SetChild(ATR, 'Property', dict(
269 DefinedOn='ClassicalEnergyHolder',
270 Name='TorsionEnergy',
271 Type='Double',
272 ))
273 SetChild(ATR, 'Property', dict(
274 DefinedOn='ClassicalEnergyHolder',
275 Name='TorsionStretchEnergy',
276 Type='Double',
277 ))
278 SetChild(ATR, 'Property', dict(
279 DefinedOn='ClassicalEnergyHolder',
280 Name='ValenceCrossTermEnergy',
281 Type='Double',
282 ))
283 SetChild(ATR, 'Property', dict(
284 DefinedOn='ClassicalEnergyHolder',
285 Name='ValenceDiagonalEnergy',
286 Type='Double',
287 ))
288 SetChild(ATR, 'Property', dict(
289 DefinedOn='ClassicalEnergyHolder',
290 Name='VanDerWaalsEnergy',
291 Type='Double',
292 ))
293 SetChild(ATR, 'Property', dict(
294 DefinedOn='SymmetrySystem',
295 Name='_Stress',
296 Type='Matrix',
297 ))
298 return ATR, XSD
301def _write_xsd_html(images, connectivity=None):
302 ATR, XSD = SetBasicChilds()
303 natoms = len(images[0])
304 atom_element = images[0].get_chemical_symbols()
305 atom_cell = images[0].get_cell()
306 atom_positions = images[0].get_positions()
307 # Set up bonds
308 bonds = []
309 if connectivity is not None:
310 for i in range(connectivity.shape[0]):
311 for j in range(i + 1, connectivity.shape[0]):
312 if connectivity[i, j]:
313 bonds.append([i, j])
314 nbonds = len(bonds)
316 # non-periodic system
317 if not images[0].pbc.all():
318 Molecule = SetChild(ATR, 'Molecule', dict(
319 ID='2',
320 NumChildren=str(natoms + nbonds),
321 Name='Lattice="1.0',
322 ))
323 # writing images[0]
324 for x in range(natoms):
325 Props = dict(
326 ID=str(x + 3),
327 Name=(atom_element[x] + str(x + 1)),
328 UserID=str(x + 1),
329 DisplayStyle=CPK_or_BnS(atom_element[x]),
330 XYZ=','.join('%1.16f' % xi for xi in atom_positions[x]),
331 Components=atom_element[x],
332 )
333 bondstr = []
334 for i, bond in enumerate(bonds):
335 if x in bond:
336 bondstr.append(str(i + 3 + natoms))
337 if bondstr:
338 Props['Connections'] = ','.join(bondstr)
339 SetChild(Molecule, 'Atom3d', Props)
340 for x in range(nbonds):
341 SetChild(Molecule, 'Bond', dict(
342 ID=str(x + 3 + natoms),
343 Connects='%i,%i' % (bonds[x][0] + 3, bonds[x][1] + 3),
344 ))
345 # periodic system
346 else:
347 atom_positions = np.dot(atom_positions, np.linalg.inv(atom_cell))
348 Props = dict(
349 ID='2',
350 Mapping='3',
351 Children=','.join(map(str, range(4, natoms + nbonds + 5))),
352 Normalized='1',
353 Name='SymmSys',
354 UserID=str(natoms + 18),
355 XYZ='0.00000000000000,0.00000000000000,0.000000000000000',
356 OverspecificationTolerance='0.05',
357 PeriodicDisplayType='Original',
358 )
359 SymmSys = SetChild(ATR, 'SymmetrySystem', Props)
361 Props = dict(
362 ID=str(natoms + nbonds + 5),
363 SymmetryDefinition=str(natoms + 4),
364 ActiveSystem='2',
365 NumFamilies='1',
366 OwnsTotalConstraintMapping='1',
367 TotalConstraintMapping='3',
368 )
369 MappngSet = SetChild(SymmSys, 'MappingSet', Props)
371 Props = dict(ID=str(natoms + nbonds + 6), NumImageMappings='0')
372 MappngFamily = SetChild(MappngSet, 'MappingFamily', Props)
374 Props = dict(
375 ID=str(natoms + len(bonds) + 7),
376 Element='1,0,0,0,0,1,0,0,0,0,1,0',
377 Constraint='1,0,0,0,0,1,0,0,0,0,1,0',
378 MappedObjects=','.join(map(str, range(4, natoms + len(bonds) + 4))),
379 DefectObjects='%i,%i' % (natoms + nbonds + 4, natoms + nbonds + 8),
380 NumImages=str(natoms + len(bonds)),
381 NumDefects='2',
382 )
383 IdentMappng = SetChild(MappngFamily, 'IdentityMapping', Props)
385 SetChild(MappngFamily, 'MappingRepairs', {'NumRepairs': '0'})
387 # writing atoms
388 for x in range(natoms):
389 Props = dict(
390 ID=str(x + 4),
391 Mapping=str(natoms + len(bonds) + 7),
392 Parent='2',
393 Name=(atom_element[x] + str(x + 1)),
394 UserID=str(x + 1),
395 DisplayStyle=CPK_or_BnS(atom_element[x]),
396 Components=atom_element[x],
397 XYZ=','.join(['%1.16f' % xi for xi in atom_positions[x]]),
398 )
399 bondstr = []
400 for i, bond in enumerate(bonds):
401 if x in bond:
402 bondstr.append(str(i + 4 * natoms + 1))
403 if bondstr:
404 Props['Connections'] = ','.join(bondstr)
405 SetChild(IdentMappng, 'Atom3d', Props)
407 for x in range(len(bonds)):
408 SetChild(IdentMappng, 'Bond', dict(
409 ID=str(x + 4 + natoms + 1),
410 Mapping=str(natoms + len(bonds) + 7),
411 Parent='2',
412 Connects='%i,%i' % (bonds[x][0] + 4, bonds[x][1] + 4),
413 ))
415 Props = dict(
416 ID=str(natoms + 4),
417 Parent='2',
418 Children=str(natoms + len(bonds) + 8),
419 DisplayStyle='Solid',
420 XYZ='0.00,0.00,0.00',
421 Color='0,0,0,0',
422 AVector=','.join(['%1.16f' % atom_cell[0, x] for x in range(3)]),
423 BVector=','.join(['%1.16f' % atom_cell[1, x] for x in range(3)]),
424 CVector=','.join(['%1.16f' % atom_cell[2, x] for x in range(3)]),
425 OrientationBase='C along Z, B in YZ plane',
426 Centering='3D Primitive-Centered',
427 Lattice='3D Triclinic',
428 GroupName='GroupName',
429 Operators='1,0,0,0,0,1,0,0,0,0,1,0',
430 DisplayRange='0,1,0,1,0,1',
431 LineThickness='2',
432 CylinderRadius='0.2',
433 LabelAxes='1',
434 ActiveSystem='2',
435 ITNumber='1',
436 LongName='P 1',
437 Qualifier='Origin-1',
438 SchoenfliesName='C1-1',
439 System='Triclinic',
440 Class='1',
441 )
442 SetChild(IdentMappng, 'SpaceGroup', Props)
444 SetChild(IdentMappng, 'ReciprocalLattice3D', dict(
445 ID=str(natoms + len(bonds) + 8),
446 Parent=str(natoms + 4),
447 ))
449 SetChild(MappngSet, 'InfiniteMapping', dict(
450 ID='3',
451 Element='1,0,0,0,0,1,0,0,0,0,1,0',
452 MappedObjects='2',
453 ))
455 return XSD, ATR
458@writer
459def write_xsd(fd, images, connectivity=None):
460 """Takes Atoms object, and write materials studio file
461 atoms: Atoms object
462 filename: path of the output file
463 connectivity: number of atoms by number of atoms matrix for connectivity
464 between atoms (0 not connected, 1 connected)
466 note: material studio file cannot use a partial periodic system. If partial
467 perodic system was inputted, full periodicity was assumed.
468 """
469 if hasattr(images, 'get_positions'):
470 images = [images]
472 XSD, ATR = _write_xsd_html(images, connectivity)
474 # Return a pretty-printed XML string for the Element.
475 rough_string = ET.tostring(XSD, 'utf-8')
476 reparsed = minidom.parseString(rough_string)
477 Document = reparsed.toprettyxml(indent='\t')
479 fd.write(Document)