Context Navigation

source: trunk/imports/G2phase_CIF.py @ 3210

Last change on this file since 3210 was 3210, checked in by vondreele, 6 years ago
some fixes for magnetic incommensurate cifs - dealing with 2 different dictionary defns for mcif!
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 42.2 KB

Line
1	# -- coding: utf-8 --
2	########### SVN repository information ###################
3	# $Date: 2017-12-29 18:52:35 +0000 (Fri, 29 Dec 2017) $
4	# $Author: vondreele $
5	# $Revision: 3210 $
6	# $URL: trunk/imports/G2phase_CIF.py $
7	# $Id: G2phase_CIF.py 3210 2017-12-29 18:52:35Z vondreele $
8	########### SVN repository information ###################
9	'''
10	Module G2phase_CIF: Coordinates from CIF
11	------------------------------------------
12
13	Parses a CIF using PyCifRW from James Hester and pulls out the
14	structural information.
15
16	If a CIF generated by ISODISTORT is encountered, extra information is
17	added to the phase entry and constraints are generated.
18
19	'''
20	# Routines to import Phase information from CIF files
21	from __future__ import division, print_function
22	import sys
23	import random as ran
24	import numpy as np
25	import re
26	import GSASIIIO as G2IO
27	import GSASIIobj as G2obj
28	import GSASIIspc as G2spc
29	import GSASIIElem as G2elem
30	import GSASIIlattice as G2lat
31	import GSASIIpy3 as G2p3
32	import GSASIIpath
33	GSASIIpath.SetVersionNumber("$Revision: 3210 $")
34	import CifFile as cif # PyCifRW from James Hester
35
36	class CIFPhaseReader(G2obj.ImportPhase):
37	'Implements a phase importer from a possibly multi-block CIF file'
38	def __init__(self):
39	super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
40	extensionlist=('.CIF','.cif','.mcif'),
41	strictExtension=False,
42	formatName = 'CIF',
43	longFormatName = 'Crystallographic Information File import'
44	)
45
46	def ContentsValidator(self, filename):
47	fp = open(filename,'r')
48	return self.CIFValidator(fp)
49	fp.close()
50
51	def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused):
52	self.isodistort_warnings = ''
53	self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
54	# make sure the ranId is really unique!
55	while self.Phase['ranId'] in usedRanIdList:
56	self.Phase['ranId'] = ran.randint(0,sys.maxsize)
57	self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
58	# make sure the ranId is really unique!
59	while self.MPhase['ranId'] in usedRanIdList:
60	self.MPhase['ranId'] = ran.randint(0,sys.maxsize)
61	returnstat = False
62	cellitems = (
63	'_cell_length_a','_cell_length_b','_cell_length_c',
64	'_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',)
65	# cellwaveitems = (
66	# '_cell_wave_vector_seq_id',
67	# '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z')
68	reqitems = (
69	'_atom_site_fract_x',
70	'_atom_site_fract_y',
71	'_atom_site_fract_z',
72	)
73	phasenamefields = (
74	'_chemical_name_common',
75	'_pd_phase_name',
76	'_chemical_formula_sum'
77	)
78	try:
79	cf = G2obj.ReadCIF(filename)
80	except cif.StarError as msg:
81	msg = 'Unreadable cif file\n'+str(msg)
82	self.errors = msg
83	self.warnings += msg
84	return False
85	# scan blocks for structural info
86	self.errors = 'Error during scan of blocks for datasets'
87	str_blklist = []
88	for blk in cf.keys():
89	for r in reqitems+cellitems:
90	if r not in cf[blk].keys():
91	break
92	else:
93	str_blklist.append(blk)
94	if not str_blklist:
95	selblk = None # no block to choose
96	elif len(str_blklist) == 1: # only one choice
97	selblk = 0
98	else: # choose from options
99	choice = []
100	for blknm in str_blklist:
101	choice.append('')
102	# accumumlate some info about this phase
103	choice[-1] += blknm + ': '
104	for i in phasenamefields: # get a name for the phase
105	name = cf[blknm].get(i,'phase name').strip()
106	if name is None or name == '?' or name == '.':
107	continue
108	else:
109	choice[-1] += name.strip()[:20] + ', '
110	break
111	na = len(cf[blknm].get("_atom_site_fract_x"))
112	if na == 1:
113	choice[-1] += '1 atom'
114	else:
115	choice[-1] += ('%d' % na) + ' atoms'
116	choice[-1] += ', cell: '
117	fmt = "%.2f,"
118	for i,key in enumerate(cellitems):
119	if i == 3: fmt = "%.f,"
120	if i == 5: fmt = "%.f"
121	choice[-1] += fmt % cif.get_number_with_esd(
122	cf[blknm].get(key))[0]
123	sg = cf[blknm].get("_symmetry_space_group_name_H-M",'')
124	if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'')
125	if not sg: sg = cf[blknm].get("_space_group_ssg_name",'')
126	if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'')
127	if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'')
128	#how about checking for super/magnetic ones as well? - reject 'X'?
129	sg = sg.replace('_','')
130	if sg: choice[-1] += ', (' + sg.strip() + ')'
131	selblk = G2IO.PhaseSelector(choice,ParentFrame=ParentFrame,
132	title= 'Select a phase from one the CIF data_ blocks below',size=(600,100))
133	self.errors = 'Error during reading of selected block'
134	#process selected phase
135	if selblk is None:
136	returnstat = False # no block selected or available
137	else: #do space group symbol & phase type first
138	blknm = str_blklist[selblk]
139	blk = cf[str_blklist[selblk]]
140	E = True
141	Super = False
142	magnetic = False
143	moddim = int(blk.get("_cell_modulation_dimension",'0'))
144	if moddim: #incommensurate
145	if moddim > 1:
146	msg = 'more than 3+1 super symmetry is not allowed in GSAS-II'
147	self.errors = msg
148	self.warnings += '\n'+msg
149	return False
150	if blk.get('_cell_subsystems_number'):
151	msg = 'Composite super structures not allowed in GSAS-II'
152	self.errors = msg
153	self.warnings += '\n'+msg
154	return False
155	sspgrp = blk.get("_space_group_ssg_name",'')
156	if not sspgrp: #might be incommensurate magnetic
157	MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'')
158	if not MSSpGrp:
159	MSSpGrp = blk.get("_space_group.magn_ssg_name",'')
160	if not MSSpGrp:
161	msg = 'No incommensurate space group name was found in the CIF.'
162	self.errors = msg
163	self.warnings += '\n'+msg
164	return False
165	if 'X' in MSSpGrp:
166	msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II'
167	self.warnings += '\n'+msg
168	self.errors = msg
169	return False
170	magnetic = True
171	if 'X' in sspgrp:
172	msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II'
173	self.warnings += '\n'+msg
174	self.errors = msg
175	return False
176	Super = True
177	if magnetic:
178	sspgrp = MSSpGrp.split('(')
179	sspgrp[1] = "("+sspgrp[1]
180	SpGrp = G2spc.StandardizeSpcName(sspgrp[0])
181	self.MPhase['General']['Type'] = 'magnetic'
182	self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
183	else:
184	sspgrp = sspgrp.split('(')
185	sspgrp[1] = "("+sspgrp[1]
186	SpGrp = sspgrp[0]
187	SpGrp = G2spc.StandardizeSpcName(SpGrp)
188	self.Phase['General']['Type'] = 'nuclear'
189	if not SpGrp:
190	print (sspgrp)
191	self.warnings += 'No space group name was found in the CIF.'
192	return False
193	SuperSg = sspgrp[1].replace('\\','')
194	SuperVec = [[0,0,.1],False,4]
195	else: #not incommensurate
196	SpGrp = blk.get("_symmetry_space_group_name_H-M",'')
197	if not SpGrp:
198	SpGrp = blk.get("_space_group_name_H-M_alt",'')
199	if not SpGrp: #try magnetic
200	MSpGrp = blk.get("_space_group.magn_name_BNS",'')
201	if not MSpGrp:
202	MSpGrp = blk.get("_space_group_magn.name_BNS",'')
203	if not MSpGrp:
204	msg = 'No recognizable space group name was found in the CIF.'
205	self.errors = msg
206	self.warnings += '\n'+msg
207	return False
208	SpGrp = MSpGrp.replace("'",'')
209	SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','') #get rid of screw '_'
210	if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:]
211	SpGrp = G2spc.StandardizeSpcName(SpGrp)
212	magnetic = True
213	self.MPhase['General']['Type'] = 'magnetic'
214	self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
215	if not SpGrp:
216	print (MSpGrp)
217	self.warnings += 'No space group name was found in the CIF.'
218	return False
219	else:
220	SpGrp = SpGrp.replace('_','')
221	self.Phase['General']['Type'] = 'nuclear'
222	#process space group symbol
223	E,SGData = G2spc.SpcGroup(SpGrp)
224	if E and SpGrp:
225	SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
226	if SpGrpNorm:
227	E,SGData = G2spc.SpcGroup(SpGrpNorm)
228	# nope, try the space group "out of the Box"
229	if E:
230	self.warnings += 'ERROR in space group symbol '+SpGrp
231	self.warnings += '\nThe space group has been set to "P 1". '
232	self.warnings += "Change this in phase's General tab."
233	self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
234	self.warnings += G2spc.SGErrors(E)
235	SGData = G2obj.P1SGData # P 1
236	self.Phase['General']['SGData'] = SGData
237
238	if magnetic and not Super:
239	SGData['SGFixed'] = True
240	try:
241	sgoploop = blk.GetLoop('_space_group_symop_magn.id')
242	sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
243	opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
244	centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]
245	except KeyError: #old mag cif names
246	sgoploop = blk.GetLoop('_space_group_symop.magn_id')
247	sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id')
248	opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1]
249	centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1]
250	SGData['SGOps'] = []
251	SGData['SGCen'] = []
252	spnflp = []
253	for op in sgoploop:
254	M,T,S = G2spc.MagText2MTS(op[opid])
255	SGData['SGOps'].append([np.array(M,dtype=float),T])
256	spnflp.append(S)
257	censpn = []
258	for cent in sgcenloop:
259	M,C,S = G2spc.MagText2MTS(cent[centid])
260	SGData['SGCen'].append(C)
261	censpn += list(np.array(spnflp)*S)
262	self.MPhase['General']['SGData'] = SGData
263	self.MPhase['General']['SGData']['SpnFlp'] = censpn
264	self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp
265	MagPtGp = blk.get('_space_group.magn_point_group')
266	if not MagPtGp:
267	MagPtGp = blk.get('_space_group_magn.point_group_name')
268	self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp
269	# GenSym,GenFlg = G2spc.GetGenSym(SGData)
270	# self.MPhase['General']['SGData']['GenSym'] = GenSym
271	# self.MPhase['General']['SGData']['GenFlg'] = GenFlg
272
273	if Super:
274	E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
275	if E:
276	self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
277	self.warnings += '\n'+E
278	SuperSg = SuperSg[:SuperSg.index(')')+1]
279	self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
280	E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
281	self.Phase['General']['SSGData'] = SSGData
282	if magnetic:
283	self.MPhase['General']['SGData'] = SGData
284	# self.MPhase['General']['SGData']['SpnFlp'] = censpn
285	self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','')
286	self.MPhase['General']['SGData']['MagPtGp'] = blk.get('_space_group.magn_point_group')
287	self.MPhase['General']['SSGData'] = SSGData
288	# GenSym,GenFlg = G2spc.GetGenSym(SGData)
289	# self.MPhase['General']['SGData']['GenSym'] = GenSym
290	# self.MPhase['General']['SGData']['GenFlg'] = GenFlg
291
292
293	# cell parameters
294	cell = []
295	for lbl in cellitems:
296	cell.append(cif.get_number_with_esd(blk[lbl])[0])
297	Volume = G2lat.calc_V(G2lat.cell2A(cell))
298	self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
299	if magnetic:
300	self.MPhase['General']['Cell'] = [False,]+cell+[Volume,]
301	if Super:
302	waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
303	waveDict = dict(waveloop.items())
304	SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
305	cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
306	cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
307
308	# read in atoms
309	self.errors = 'Error during reading of atoms'
310	atomlbllist = [] # table to look up atom IDs
311	atomloop = blk.GetLoop('_atom_site_label')
312	atomkeys = [i.lower() for i in atomloop.keys()]
313	if not blk.get('_atom_site_type_symbol'):
314	self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
315	if magnetic:
316	try:
317	magmoment = '_atom_site_moment.label'
318	magatomloop = blk.GetLoop(magmoment)
319	magatomkeys = [i.lower() for i in magatomloop.keys()]
320	magatomlabels = blk.get(magmoment)
321	G2MagDict = {'_atom_site_moment.label': 0,
322	'_atom_site_moment.crystalaxis_x':7,
323	'_atom_site_moment.crystalaxis_y':8,
324	'_atom_site_moment.crystalaxis_z':9}
325	except KeyError:
326	magmoment = '_atom_site_moment_label'
327	magatomloop = blk.GetLoop(magmoment)
328	magatomkeys = [i.lower() for i in magatomloop.keys()]
329	magatomlabels = blk.get(magmoment)
330	G2MagDict = {'_atom_site_moment_label': 0,
331	'_atom_site_moment_crystalaxis_x':7,
332	'_atom_site_moment_crystalaxis_y':8,
333	'_atom_site_moment_crystalaxis_z':9}
334
335	if blk.get('_atom_site_aniso_label'):
336	anisoloop = blk.GetLoop('_atom_site_aniso_label')
337	anisokeys = [i.lower() for i in anisoloop.keys()]
338	anisolabels = blk.get('_atom_site_aniso_label')
339	else:
340	anisoloop = None
341	anisokeys = []
342	anisolabels = []
343	if Super:
344	occFloop = None
345	occCloop = None
346	occFdict = {}
347	occCdict = {}
348	displSloop = None
349	displFloop = None
350	MagFloop = None
351	displSdict = {}
352	displFdict = {}
353	MagFdict = {}
354	UijFloop = None
355	UijFdict = {}
356	if blk.get('_atom_site_occ_Fourier_atom_site_label'):
357	occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
358	occFdict = dict(occFloop.items())
359	if blk.get('_atom_site_occ_special_func_atom_site_label'): #Crenel (i.e. Block Wave) occ
360	occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
361	occCdict = dict(occCloop.items())
362	if blk.get('_atom_site_displace_Fourier_atom_site_label'):
363	displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
364	displFdict = dict(displFloop.items())
365	if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
366	displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
367	displSdict = dict(displSloop.items())
368	if blk.get('_atom_site_U_Fourier_atom_site_label'):
369	UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
370	UijFdict = dict(UijFloop.items())
371	if blk.get('_atom_site_moment_Fourier_atom_site_label'):
372	MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label')
373	MagFdict = dict(MagFloop.items())
374	Mnames = ['_atom_site_moment_fourier_atom_site_label',
375	'_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id',
376	'_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos']
377	elif blk.get('_atom_site_moment_Fourier.atom_site_label'):
378	MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label')
379	MagFdict = dict(MagFloop.items())
380	Mnames = ['_atom_site_moment_fourier.atom_site_label',
381	'_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id',
382	'_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos']
383	self.Phase['Atoms'] = []
384	if magnetic:
385	self.MPhase['Atoms'] = []
386	G2AtomDict = { '_atom_site_type_symbol' : 1,
387	'_atom_site_label' : 0,
388	'_atom_site_fract_x' : 3,
389	'_atom_site_fract_y' : 4,
390	'_atom_site_fract_z' : 5,
391	'_atom_site_occupancy' : 6,
392	'_atom_site_aniso_u_11' : 11,
393	'_atom_site_aniso_u_22' : 12,
394	'_atom_site_aniso_u_33' : 13,
395	'_atom_site_aniso_u_12' : 14,
396	'_atom_site_aniso_u_13' : 15,
397	'_atom_site_aniso_u_23' : 16, }
398
399	ranIdlookup = {}
400	for aitem in atomloop:
401	atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
402	for val,key in zip(aitem,atomkeys):
403	col = G2AtomDict.get(key,-1)
404	if col >= 3:
405	atomlist[col] = cif.get_number_with_esd(val)[0]
406	if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
407	elif col is not None and col != -1:
408	atomlist[col] = val
409	elif key in ('_atom_site_thermal_displace_type',
410	'_atom_site_adp_type'): #Iso or Aniso?
411	if val.lower() == 'uani':
412	atomlist[9] = 'A'
413	elif key == '_atom_site_u_iso_or_equiv':
414	uisoval = cif.get_number_with_esd(val)[0]
415	if uisoval is not None:
416	atomlist[10] = uisoval
417	if not atomlist[1] and atomlist[0]:
418	typ = atomlist[0].rstrip('0123456789-+')
419	if G2elem.CheckElement(typ):
420	atomlist[1] = typ
421	if not atomlist[1]:
422	atomlist[1] = 'Xe'
423	self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
424	if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
425	atomlist[9] = 'A'
426	for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
427	col = G2AtomDict.get(key)
428	if col:
429	atomlist[col] = cif.get_number_with_esd(val)[0]
430	atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
431	atomlist[1] = G2elem.FixValence(atomlist[1])
432	atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
433	self.Phase['Atoms'].append(atomlist)
434	ranIdlookup[atomlist[0]] = atomlist[-1]
435	if atomlist[0] in atomlbllist:
436	self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
437	else:
438	atomlbllist.append(atomlist[0])
439
440	if magnetic and atomlist[0] in magatomlabels:
441	matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
442	for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
443	mcol = G2MagDict.get(mkey,-1)
444	if mcol:
445	matomlist[mcol] = cif.get_number_with_esd(mval)[0]
446	self.MPhase['Atoms'].append(matomlist)
447	if Super:
448	Sfrac = []
449	Sadp = np.zeros((4,12))
450	Spos = np.zeros((4,6))
451	Smag = np.zeros((4,6))
452	nim = -1
453	waveType = 'Fourier'
454	if displFdict:
455	for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
456	if item == atomlist[0]:
457	waveType = 'Fourier'
458	ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
459	im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
460	if im != nim:
461	nim = im
462	val = displFdict['_atom_site_displace_fourier_param_sin'][i]
463	Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
464	val = displFdict['_atom_site_displace_fourier_param_cos'][i]
465	Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
466	if nim >= 0:
467	Spos = [[spos,False] for spos in Spos[:nim]]
468	else:
469	Spos = []
470	if UijFdict:
471	nim = -1
472	for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
473	if item == atomlist[0]:
474	ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
475	im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
476	if im != nim:
477	nim = im
478	val = UijFdict['_atom_site_u_fourier_param_sin'][i]
479	Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
480	val = UijFdict['_atom_site_u_fourier_param_cos'][i]
481	Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
482	if nim >= 0:
483	Sadp = [[sadp,False] for sadp in Sadp[:nim]]
484	else:
485	Sadp = []
486	if MagFdict:
487	nim = -1
488	for i,item in enumerate(MagFdict[Mnames[0]]):
489	if item == atomlist[0]:
490	waveType = 'Fourier'
491	ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
492	im = int(MagFdict[Mnames[2]][i])
493	if im != nim:
494	nim = im
495	val = MagFdict[Mnames[3]][i]
496	Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
497	val = MagFdict[Mnames[4]][i]
498	Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
499	if nim >= 0:
500	Smag = [[smag,False] for smag in Smag[:nim]]
501	else:
502	Smag = []
503	SSdict = {'SS1':{'waveType':waveType,'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
504	if magnetic and atomlist[0] in magatomlabels:
505	matomlist.append(SSdict)
506	atomlist.append(SSdict)
507	if len(atomlbllist) != len(self.Phase['Atoms']):
508	self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
509	for lbl in phasenamefields: # get a name for the phase
510	name = blk.get(lbl)
511	if name is None:
512	continue
513	name = name.strip()
514	if name == '?' or name == '.':
515	continue
516	else:
517	break
518	else: # no name found, use block name for lack of a better choice
519	name = blknm
520	self.Phase['General']['Name'] = name.strip()[:20]
521	self.Phase['General']['Super'] = Super
522	if magnetic:
523	self.MPhase['General']['Type'] = 'magnetic'
524	self.MPhase['General']['Name'] = name.strip()[:20]+' mag'
525	self.MPhase['General']['Super'] = Super
526	if Super:
527	self.MPhase['General']['Modulated'] = True
528	self.MPhase['General']['SuperVec'] = SuperVec
529	self.MPhase['General']['SuperSg'] = SuperSg
530	self.MPhase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
531	else:
532	self.MPhase = None
533	if Super:
534	self.Phase['General']['Modulated'] = True
535	self.Phase['General']['SuperVec'] = SuperVec
536	self.Phase['General']['SuperSg'] = SuperSg
537	self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
538	if not self.isodistort_warnings:
539	if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'):
540	self.errors = "Error while processing ISODISTORT constraints"
541	self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
542	else:
543	self.warnings += self.isodistort_warnings
544	returnstat = True
545	return returnstat
546
547	def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
548	'Process ISODISTORT items to create constraints etc.'
549	varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
550	'Maps ISODISTORT parm names to GSAS-II names'
551	#----------------------------------------------------------------------
552	# read in the ISODISTORT displacement modes
553	#----------------------------------------------------------------------
554	self.Constraints = []
555	explaination = {}
556	if blk.get('_iso_displacivemode_label'):
557	modelist = []
558	shortmodelist = []
559	for lbl in blk.get('_iso_displacivemode_label'):
560	modelist.append(lbl)
561	# assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
562	# where SSSSS is the parent spacegroup, [x,y,z] is a location
563	regexp = re.match(r'.?\[.?\](.?)\(.?\)(.*)',lbl)
564	# this extracts the AAAAA and BBBBB parts of the string
565	if regexp:
566	lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
567	G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
568	# read in the coordinate offset variables names and map them to G2 names/objects
569	coordVarLbl = []
570	G2varLbl = []
571	G2varObj = []
572	error = False
573	for lbl in blk.get('_iso_deltacoordinate_label'):
574	coordVarLbl.append(lbl)
575	if '_' in lbl:
576	albl = lbl[:lbl.rfind('_')]
577	vlbl = lbl[lbl.rfind('_')+1:]
578	else:
579	self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
580	error = True
581	continue
582	if albl not in atomlbllist:
583	self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
584	error = True
585	continue
586	else:
587	anum = atomlbllist.index(albl)
588	var = varLookup.get(vlbl)
589	if not var:
590	self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
591	error = True
592	continue
593	G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
594	G2varObj.append(G2obj.G2VarObj(
595	(self.Phase['ranId'],None,var,ranIdlookup[albl])
596	))
597	if error:
598	raise Exception("Error decoding variable labels")
599
600	if len(G2varObj) != len(modelist):
601	print ("non-square input")
602	raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
603
604	error = False
605	ParentCoordinates = {}
606	for lbl,exp in zip(
607	blk.get('_iso_coordinate_label'),
608	blk.get('_iso_coordinate_formula'),
609	):
610	if '_' in lbl:
611	albl = lbl[:lbl.rfind('_')]
612	vlbl = lbl[lbl.rfind('_')+1:]
613	else:
614	self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
615	error = True
616	continue
617	if vlbl not in 'xyz' or len(vlbl) != 1:
618	self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
619	error = True
620	continue
621	i = 'xyz'.index(vlbl)
622	if not ParentCoordinates.get(albl):
623	ParentCoordinates[albl] = [None,None,None]
624	if '+' in exp:
625	val = exp.split('+')[0].strip()
626	val = G2p3.FormulaEval(val)
627	if val is None:
628	self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
629	error = True
630	continue
631	ParentCoordinates[albl][i] = val
632	else:
633	ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
634	if error:
635	print (self.warnings)
636	raise Exception("Error decoding variable labels")
637	# get mapping of modes to atomic coordinate displacements
638	displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
639	for row,col,val in zip(
640	blk.get('_iso_displacivemodematrix_row'),
641	blk.get('_iso_displacivemodematrix_col'),
642	blk.get('_iso_displacivemodematrix_value'),):
643	displacivemodematrix[int(row)-1,int(col)-1] = float(val)
644	# Invert to get mapping of atom displacements to modes
645	displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
646	# create the constraints
647	for i,row in enumerate(displacivemodeInvmatrix):
648	constraint = []
649	for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
650	if k == 0: continue
651	constraint.append([k,G2varObj[j]])
652	constraint += [shortmodelist[i],False,'f']
653	self.Constraints.append(constraint)
654	#----------------------------------------------------------------------
655	# save the ISODISTORT info for "mode analysis"
656	if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
657	self.Phase['ISODISTORT'].update({
658	'IsoModeList' : modelist,
659	'G2ModeList' : shortmodelist,
660	'IsoVarList' : coordVarLbl,
661	'G2VarList' : G2varObj,
662	'ParentStructure' : ParentCoordinates,
663	'Var2ModeMatrix' : displacivemodeInvmatrix,
664	'Mode2VarMatrix' : displacivemodematrix,
665	})
666	# make explaination dictionary
667	for mode,shortmode in zip(modelist,shortmodelist):
668	explaination[shortmode] = "ISODISTORT full name "+str(mode)
669	#----------------------------------------------------------------------
670	# now read in the ISODISTORT occupancy modes
671	#----------------------------------------------------------------------
672	if blk.get('_iso_occupancymode_label'):
673	modelist = []
674	shortmodelist = []
675	for lbl in blk.get('_iso_occupancymode_label'):
676	modelist.append(lbl)
677	# assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
678	# where SSSSS is the parent spacegroup, [x,y,z] is a location
679	regexp = re.match(r'.?\[.?\](.?)\(.?\)(.*)',lbl)
680	# this extracts the AAAAA and BBBBB parts of the string
681	if regexp:
682	lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
683	lbl = lbl.replace('order','o')
684	G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
685	# read in the coordinate offset variables names and map them to G2 names/objects
686	occVarLbl = []
687	G2varLbl = []
688	G2varObj = []
689	error = False
690	for lbl in blk.get('_iso_deltaoccupancy_label'):
691	occVarLbl.append(lbl)
692	if '_' in lbl:
693	albl = lbl[:lbl.rfind('_')]
694	vlbl = lbl[lbl.rfind('_')+1:]
695	else:
696	self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
697	error = True
698	continue
699	if albl not in atomlbllist:
700	self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
701	error = True
702	continue
703	else:
704	anum = atomlbllist.index(albl)
705	var = varLookup.get(vlbl)
706	if not var:
707	self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
708	error = True
709	continue
710	G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID
711	G2varObj.append(G2obj.G2VarObj(
712	(self.Phase['ranId'],None,var,ranIdlookup[albl])
713	))
714	if error:
715	raise Exception("Error decoding variable labels")
716
717	if len(G2varObj) != len(modelist):
718	print ("non-square input")
719	raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
720
721	error = False
722	ParentCoordinates = {}
723	for lbl,exp in zip(
724	blk.get('_iso_occupancy_label'),
725	blk.get('_iso_occupancy_formula'),
726	):
727	if '_' in lbl:
728	albl = lbl[:lbl.rfind('_')]
729	vlbl = lbl[lbl.rfind('_')+1:]
730	else:
731	self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
732	error = True
733	continue
734	if vlbl != 'occ':
735	self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
736	error = True
737	continue
738	if '+' in exp:
739	val = exp.split('+')[0].strip()
740	val = G2p3.FormulaEval(val)
741	if val is None:
742	self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
743	error = True
744	continue
745	ParentCoordinates[albl] = val
746	if error:
747	raise Exception("Error decoding occupancy labels")
748	# get mapping of modes to atomic coordinate displacements
749	occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
750	for row,col,val in zip(
751	blk.get('_iso_occupancymodematrix_row'),
752	blk.get('_iso_occupancymodematrix_col'),
753	blk.get('_iso_occupancymodematrix_value'),):
754	occupancymodematrix[int(row)-1,int(col)-1] = float(val)
755	# Invert to get mapping of atom displacements to modes
756	occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
757	# create the constraints
758	for i,row in enumerate(occupancymodeInvmatrix):
759	constraint = []
760	for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
761	if k == 0: continue
762	constraint.append([k,G2varObj[j]])
763	constraint += [shortmodelist[i],False,'f']
764	self.Constraints.append(constraint)
765	#----------------------------------------------------------------------
766	# save the ISODISTORT info for "mode analysis"
767	if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
768	self.Phase['ISODISTORT'].update({
769	'OccModeList' : modelist,
770	'G2OccModeList' : shortmodelist,
771	'OccVarList' : occVarLbl,
772	'G2OccVarList' : G2varObj,
773	'BaseOcc' : ParentCoordinates,
774	'Var2OccMatrix' : occupancymodeInvmatrix,
775	'Occ2VarMatrix' : occupancymodematrix,
776	})
777	# make explaination dictionary
778	for mode,shortmode in zip(modelist,shortmodelist):
779	explaination[shortmode] = "ISODISTORT full name "+str(mode)
780	#----------------------------------------------------------------------
781	# done with read
782	#----------------------------------------------------------------------
783	if explaination: self.Constraints.append(explaination)
784
785	# # debug: show the mode var to mode relations
786	# for i,row in enumerate(displacivemodeInvmatrix):
787	# l = ''
788	# for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
789	# if k == 0: continue
790	# if l: l += ' + '
791	# #l += lbl+' * '+str(k)
792	# l += G2varLbl[j]+' * '+str(k)
793	# print str(i) + ': '+shortmodelist[i]+' = '+l
794	# print 70*'='
795
796	# # debug: Get the ISODISTORT offset values
797	# coordVarDelta = {}
798	# for lbl,val in zip(
799	# blk.get('_iso_deltacoordinate_label'),
800	# blk.get('_iso_deltacoordinate_value'),):
801	# coordVarDelta[lbl] = float(val)
802	# modeVarDelta = {}
803	# for lbl,val in zip(
804	# blk.get('_iso_displacivemode_label'),
805	# blk.get('_iso_displacivemode_value'),):
806	# modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
807
808	# print 70*'='
809	# # compute the mode values from the reported coordinate deltas
810	# for i,row in enumerate(displacivemodeInvmatrix):
811	# l = ''
812	# sl = ''
813	# s = 0.
814	# for lbl,k in zip(coordVarLbl,row):
815	# if k == 0: continue
816	# if l: l += ' + '
817	# l += lbl+' * '+str(k)
818	# if sl: sl += ' + '
819	# sl += str(coordVarDelta[lbl])+' * '+str(k)
820	# s += coordVarDelta[lbl] * k
821	# print 'a'+str(i)+' = '+l
822	# print '\t= '+sl
823	# print modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s
824	# print
825
826	# print 70*'='
827	# # compute the coordinate displacements from the reported mode values
828	# for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
829	# l = ''
830	# sl = ''
831	# s = 0.0
832	# for j,k in enumerate(row):
833	# if k == 0: continue
834	# if l: l += ' + '
835	# l += 'a'+str(j+1)+' * '+str(k)
836	# if sl: sl += ' + '
837	# sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k)
838	# s += modeVarDelta[modelist[j]] * k
839	# print lbl+' = '+l
840	# print '\t= '+sl
841	# print lbl,G2varLbl[i],coordVarDelta[lbl],s
842	# print
843
844	# determine the coordinate delta values from deviations from the parent structure
845	# for atmline in self.Phase['Atoms']:
846	# lbl = atmline[0]
847	# x,y,z = atmline[3:6]
848	# if lbl not in ParentCoordinates:
849	# print lbl,x,y,z
850	# continue
851	# px,py,pz = ParentCoordinates[lbl]
852	# print lbl,x,y,z,x-px,y-py,z-pz

Note: See TracBrowser for help on using the repository browser.

Download in other formats: