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