source: trunk/imports/G2phase_CIF.py @ 4475

Last change on this file since 4475 was 4475, checked in by vondreele, 3 years ago

fix prolem in cif files where anisotropic params are as ...... --> None; atom is now set to 'I' instead of 'A'
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 58.7 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2020-06-09 15:36:09 +0000 (Tue, 09 Jun 2020) $
4# $Author: vondreele $
5# $Revision: 4475 $
6# $URL: trunk/imports/G2phase_CIF.py $
7# $Id: G2phase_CIF.py 4475 2020-06-09 15:36:09Z vondreele $
8########### SVN repository information ###################
9'''
10*Module G2phase_CIF: Coordinates from CIF*
11------------------------------------------
12
13Parses a CIF using  PyCifRW from James Hester and pulls out the
14structural information.
15
16If a CIF generated by ISODISTORT is encountered, extra information is
17added to the phase entry and constraints are generated.
18
19'''
20# Routines to import Phase information from CIF files
21from __future__ import division, print_function
22import sys
23import random as ran
24import numpy as np
25import re
26import copy
27import GSASIIobj as G2obj
28import GSASIIspc as G2spc
29import GSASIIElem as G2elem
30import GSASIIlattice as G2lat
31import GSASIIpy3 as G2p3
32import GSASIIpath
33try:
34    import GSASIIctrlGUI as G2G
35except ImportError:
36    pass
37GSASIIpath.SetVersionNumber("$Revision: 4475 $")
38import CifFile as cif # PyCifRW from James Hester
39debug = GSASIIpath.GetConfigValue('debug')
40#debug = False
41
42class CIFPhaseReader(G2obj.ImportPhase):
43    'Implements a phase importer from a possibly multi-block CIF file'
44    def __init__(self):
45        super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__
46            extensionlist=('.CIF','.cif','.mcif'),
47            strictExtension=False,
48            formatName = 'CIF',
49            longFormatName = 'Crystallographic Information File import'
50            )
51       
52    def ContentsValidator(self, filename):
53        fp = open(filename,'r')
54        return self.CIFValidator(fp)
55        fp.close()
56
57    def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused):
58        isodistort_warnings = '' # errors that would prevent an isodistort analysis
59        self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
60        # make sure the ranId is really unique!
61        while self.Phase['ranId'] in usedRanIdList:
62            self.Phase['ranId'] = ran.randint(0,sys.maxsize)
63        self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict
64        # make sure the ranId is really unique!
65        while self.MPhase['ranId'] in usedRanIdList:
66            self.MPhase['ranId'] = ran.randint(0,sys.maxsize)
67        returnstat = False
68        cellitems = (
69            '_cell_length_a','_cell_length_b','_cell_length_c',
70            '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',)
71#        cellwaveitems = (
72#            '_cell_wave_vector_seq_id',
73#            '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z')
74        reqitems = (
75             '_atom_site_fract_x',
76             '_atom_site_fract_y',
77             '_atom_site_fract_z',
78            )
79        phasenamefields = (
80            '_chemical_name_common',
81            '_pd_phase_name',
82            '_chemical_formula_sum'
83            )
84        try:
85            cf = G2obj.ReadCIF(filename)
86        except cif.StarError as msg:
87            msg  = 'Unreadable cif file\n'+str(msg)
88            self.errors = msg
89            self.warnings += msg
90            return False
91        # scan blocks for structural info
92        self.errors = 'Error during scan of blocks for datasets'
93        str_blklist = []
94        for blk in cf.keys():
95            for r in reqitems+cellitems:
96                if r not in cf[blk].keys():
97                    break
98            else:
99                str_blklist.append(blk)
100        if not str_blklist:
101            selblk = None # no block to choose
102        elif len(str_blklist) == 1: # only one choice
103            selblk = 0
104        else:                       # choose from options
105            choice = []
106            for blknm in str_blklist:
107                choice.append('')
108                # accumumlate some info about this phase
109                choice[-1] += blknm + ': '
110                for i in phasenamefields: # get a name for the phase
111                    name = cf[blknm].get(i,'phase name').strip()
112                    if name is None or name == '?' or name == '.':
113                        continue
114                    else:
115                        choice[-1] += name.strip() + ', '
116                        break
117                na = len(cf[blknm].get("_atom_site_fract_x"))
118                if na == 1:
119                    choice[-1] += '1 atom'
120                else:
121                    choice[-1] += ('%d' % na) + ' atoms'
122                choice[-1] += ', cell: '
123                fmt = "%.2f,"
124                for i,key in enumerate(cellitems):
125                    if i == 3: fmt = "%.f,"
126                    if i == 5: fmt = "%.f"
127                    choice[-1] += fmt % cif.get_number_with_esd(
128                        cf[blknm].get(key))[0]
129                sg = cf[blknm].get("_symmetry_space_group_name_H-M",'')
130                if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'')
131                if not sg: sg = cf[blknm].get("_space_group_ssg_name",'')
132                if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'')
133                if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'')
134                #how about checking for super/magnetic ones as well? - reject 'X'?
135                sg = sg.replace('_','')
136                if sg: choice[-1] += ', (' + sg.strip() + ')'
137            selblk = G2G.PhaseSelector(choice,ParentFrame=ParentFrame,
138                title= 'Select a phase from one the CIF data_ blocks below',size=(600,100))
139        self.errors = 'Error during reading of selected block'
140#process selected phase
141        if selblk is None:
142            returnstat = False # no block selected or available
143        else:   #do space group symbol & phase type first
144            blknm = str_blklist[selblk]
145            blk = cf[str_blklist[selblk]]
146            E = True
147            Super = False
148            magnetic = False
149            moddim = int(blk.get("_cell_modulation_dimension",'0'))
150            if moddim:      #incommensurate
151                if moddim > 1:
152                    msg = 'more than 3+1 super symmetry is not allowed in GSAS-II'
153                    self.errors = msg
154                    self.warnings += '\n'+msg
155                    return False
156                if blk.get('_cell_subsystems_number'):
157                    msg = 'Composite super structures not allowed in GSAS-II'
158                    self.errors = msg
159                    self.warnings += '\n'+msg
160                    return False
161                sspgrp = blk.get("_space_group_ssg_name",'')
162                if not sspgrp:          #might be incommensurate magnetic
163                    MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'')
164                    if not MSSpGrp:
165                        MSSpGrp = blk.get("_space_group.magn_ssg_name",'') 
166                    if not MSSpGrp:
167                        msg = 'No incommensurate space group name was found in the CIF.'
168                        self.errors = msg
169                        self.warnings += '\n'+msg
170                        return False                                                           
171                    if 'X' in MSSpGrp:
172                        msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II'
173                        self.warnings += '\n'+msg
174                        self.errors = msg
175                        return False
176                    magnetic = True
177                if 'X' in sspgrp:
178                    msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II'
179                    self.warnings += '\n'+msg
180                    self.errors = msg
181                    return False
182                Super = True
183                if magnetic:
184                    sspgrp = MSSpGrp.split('(')
185                    sspgrp[1] = "("+sspgrp[1]
186                    SpGrp = G2spc.StandardizeSpcName(sspgrp[0])
187                    if "1'" in SpGrp: sspgrp[1] = sspgrp[1][:-1]  #take off extra 's'; gets put back later
188                    MSpGrp = sspgrp[0]
189                    self.MPhase['General']['Type'] = 'magnetic'
190                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
191                else:
192                    sspgrp = sspgrp.split('(')
193                    sspgrp[1] = "("+sspgrp[1]
194                    SpGrp = sspgrp[0]
195                    if "1'" in SpGrp:       #nonmagnetics can't be gray
196                        SpGrp = SpGrp.replace("1'",'')
197                        sspgrp[1] = sspgrp[1][:-1]  #take off extra 's'
198#                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
199                    self.Phase['General']['Type'] = 'nuclear'
200                if not SpGrp:
201                    print (sspgrp)
202                    self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II'
203                    return False
204                SuperSg = sspgrp[1].replace('\\','').replace(',','')
205                SuperVec = [[0,0,.1],False,4]
206            else:   #not incommensurate
207                SpGrp = blk.get("_symmetry_space_group_name_H-M",'')
208                if not SpGrp:
209                    SpGrp = blk.get("_space_group_name_H-M_alt",'')
210                if not SpGrp:   #try magnetic           
211                    MSpGrp = blk.get("_space_group.magn_name_BNS",'')
212                    if not MSpGrp:
213                        MSpGrp = blk.get("_space_group_magn.name_BNS",'')
214                        if not MSpGrp:
215                            msg = 'No recognizable space group name was found in the CIF.'
216                            self.errors = msg
217                            self.warnings += '\n'+msg
218                            return False
219                    SpGrp = blk.get('_parent_space_group.name_H-M_alt')
220                    if not SpGrp:
221                        SpGrp = blk.get('_parent_space_group.name_H-M')
222#                    SpGrp = MSpGrp.replace("'",'')
223                    SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','')   #get rid of screw '_'
224                    if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:]
225                    SpGrp = G2spc.StandardizeSpcName(SpGrp)
226                    magnetic = True
227                    self.MPhase['General']['Type'] = 'magnetic'
228                    self.MPhase['General']['AtomPtrs'] = [3,1,10,12]
229                    if not SpGrp:
230                        print (MSpGrp)
231                        self.warnings += 'No space group name was found in the CIF.'
232                        return False
233                else:
234                    SpGrp = SpGrp.replace('_','')
235                    self.Phase['General']['Type'] = 'nuclear'
236#process space group symbol
237            E,SGData = G2spc.SpcGroup(SpGrp)
238            if E and SpGrp:
239                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
240                if SpGrpNorm:
241                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
242            # nope, try the space group "out of the Box"
243            if E:
244                self.warnings += 'ERROR in space group symbol '+SpGrp
245                self.warnings += '\nThe space group has been set to "P 1". '
246                self.warnings += "Change this in phase's General tab."
247                self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: '
248                self.warnings += G2spc.SGErrors(E)
249                SGData = G2obj.P1SGData # P 1
250            self.Phase['General']['SGData'] = SGData
251            if Super:
252                E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
253                if E:
254                    self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg
255                    self.warnings += '\n'+E
256                    SuperSg = SuperSg[:SuperSg.index(')')+1]
257                    self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg
258                    E,SSGData = G2spc.SSpcGroup(SGData,SuperSg)
259                self.Phase['General']['SSGData'] = SSGData
260                if magnetic:
261                    self.MPhase['General']['SGData'] = SGData
262                    self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','')
263                    self.MPhase['General']['SSGData'] = SSGData
264
265            if magnetic:    #replace std operaors with those from cif file - probably not the same!
266                SGData['SGFixed'] = True
267                SGData['SGOps'] = []
268                SGData['SGCen'] = []
269                if Super:
270                    SSGData['SSGOps'] = []
271                    SSGData['SSGCen'] = []
272                    try:
273                        sgoploop = blk.GetLoop('_space_group_symop_magn_ssg_operation.id')
274                        sgcenloop = blk.GetLoop('_space_group_symop_magn_ssg_centering.id')
275                        opid = sgoploop.GetItemPosition('_space_group_symop_magn_ssg_operation.algebraic')[1]
276                        centid = sgcenloop.GetItemPosition('_space_group_symop_magn_ssg_centering.algebraic')[1]                   
277                    except KeyError:        #old mag cif names
278                        sgoploop = blk.GetLoop('_space_group_symop.magn_ssg_id')
279                        sgcenloop = blk.GetLoop('_space_group_symop.magn_ssg_centering_id')
280                        opid = sgoploop.GetItemPosition('_space_group_symop.magn_ssg_operation_algebraic')[1]
281                        centid = sgcenloop.GetItemPosition('_space_group_symop.magn_ssg_centering_algebraic')[1]
282                    spnflp = []
283                    for op in sgoploop:
284                        M,T,S = G2spc.MagSSText2MTS(op[opid])
285                        SSGData['SSGOps'].append([np.array(M,dtype=float),T])
286                        SGData['SGOps'].append([np.array(M,dtype=float)[:3,:3],T[:3]])
287                        spnflp.append(S)
288                    censpn = []
289                    for cent in sgcenloop:
290                        M,C,S = G2spc.MagSSText2MTS(cent[centid])
291                        SSGData['SSGCen'].append(C)
292                        SGData['SGCen'].append(C[:3])
293                        censpn += list(np.array(spnflp)*S)
294                    self.MPhase['General']['SSGData'] = SSGData
295                else:
296                    try:
297                        sgoploop = blk.GetLoop('_space_group_symop_magn_operation.id')
298                        opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
299                        try:
300                            sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
301                            centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]
302                        except KeyError:
303                            sgcenloop = None
304                    except KeyError:                                         
305                        try:
306                            sgoploop = blk.GetLoop('_space_group_symop_magn.id')
307                            sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id')
308                            opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1]
309                            centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1]                   
310                        except KeyError:        #old mag cif names
311                            sgoploop = blk.GetLoop('_space_group_symop.magn_id')
312                            sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id')
313                            opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1]
314                            centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1]
315                    spnflp = []
316                    for op in sgoploop:
317                        try:
318                            M,T,S = G2spc.MagText2MTS(op[opid])
319                            SGData['SGOps'].append([np.array(M,dtype=float),T])
320                            spnflp.append(S)
321                        except KeyError:
322                            self.warnings += 'Space group operator '+op[opid]+' is not recognized by GSAS-II'
323                            return False
324                    censpn = []
325                    if sgcenloop:
326                        for cent in sgcenloop:
327                            M,C,S = G2spc.MagText2MTS(cent[centid])
328                            SGData['SGCen'].append(C)
329                            censpn += list(np.array(spnflp)*S)
330                    else:
331                            M,C,S = G2spc.MagText2MTS('x,y,z,+1')
332                            SGData['SGCen'].append(C)
333                            censpn += list(np.array(spnflp)*S)                       
334                self.MPhase['General']['SGData'] = SGData
335                self.MPhase['General']['SGData']['SpnFlp'] = censpn
336                G2spc.GenMagOps(SGData)         #set magMom
337                self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp
338                if "1'" in MSpGrp:
339                    SGData['SGGray'] = True
340                MagPtGp = blk.get('_space_group.magn_point_group')
341                if not MagPtGp:
342                    MagPtGp = blk.get('_space_group.magn_point_group_name')
343                if not MagPtGp:
344                    MagPtGp = blk.get('_space_group_magn.point_group_name')
345                self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp
346
347            # cell parameters
348            cell = []
349            for lbl in cellitems:
350                cell.append(cif.get_number_with_esd(blk[lbl])[0])
351            Volume = G2lat.calc_V(G2lat.cell2A(cell))
352            self.Phase['General']['Cell'] = [False,]+cell+[Volume,]
353            if magnetic:
354                self.MPhase['General']['Cell'] = [False,]+cell+[Volume,]               
355            if Super:
356                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
357                waveDict = dict(waveloop.items())
358                SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
359                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
360                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
361
362            # read in atoms
363            self.errors = 'Error during reading of atoms'
364            atomlbllist = [] # table to look up atom IDs
365            atomloop = blk.GetLoop('_atom_site_label')
366            atomkeys = [i.lower() for i in atomloop.keys()]
367            if not blk.get('_atom_site_type_symbol'):
368                isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed'
369            if magnetic:
370                try:
371                    magmoment = '_atom_site_moment.label'
372                    magatomloop = blk.GetLoop(magmoment)
373                    magatomkeys = [i.lower() for i in magatomloop.keys()]
374                    magatomlabels = blk.get(magmoment)
375                    G2MagDict = {'_atom_site_moment.label': 0,
376                                 '_atom_site_moment.crystalaxis_x':7,
377                                 '_atom_site_moment.crystalaxis_y':8,
378                                 '_atom_site_moment.crystalaxis_z':9}
379                except KeyError:
380                    magmoment = '_atom_site_moment_label'
381                    magatomloop = blk.GetLoop(magmoment)
382                    magatomkeys = [i.lower() for i in magatomloop.keys()]
383                    magatomlabels = blk.get(magmoment)
384                    G2MagDict = {'_atom_site_moment_label': 0,
385                                 '_atom_site_moment_crystalaxis_x':7,
386                                 '_atom_site_moment_crystalaxis_y':8,
387                                 '_atom_site_moment_crystalaxis_z':9}
388                   
389            if blk.get('_atom_site_aniso_label'):
390                anisoloop = blk.GetLoop('_atom_site_aniso_label')
391                anisokeys = [i.lower() for i in anisoloop.keys()]
392                anisolabels = blk.get('_atom_site_aniso_label')
393            else:
394                anisoloop = None
395                anisokeys = []
396                anisolabels = []
397            if Super:
398                occFloop = None
399                occCloop = None
400                occFdict = {}
401                occCdict = {}
402                displSloop = None
403                displFloop = None
404                MagFloop = None
405                displSdict = {}
406                displFdict = {}
407                MagFdict = {}
408                UijFloop = None
409                UijFdict = {}
410                #occupancy modulation
411                if blk.get('_atom_site_occ_Fourier_atom_site_label'):
412                    occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label')
413                    occFdict = dict(occFloop.items())
414                if blk.get('_atom_site_occ_special_func_atom_site_label'):  #Crenel (i.e. Block Wave) occ
415                    occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label')
416                    occCdict = dict(occCloop.items())
417                #position modulation
418                if blk.get('_atom_site_displace_Fourier_atom_site_label'):
419                    displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label')
420                    displFdict = dict(displFloop.items())                           
421                if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth
422                    displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label')
423                    displSdict = dict(displSloop.items())
424                #U modulation
425                if blk.get('_atom_site_U_Fourier_atom_site_label'):
426                    UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label')
427                    UijFdict = dict(UijFloop.items())
428                #Mag moment modulation
429                if blk.get('_atom_site_moment_Fourier_atom_site_label'):
430                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label')
431                    MagFdict = dict(MagFloop.items())
432                    Mnames =  ['_atom_site_moment_fourier_atom_site_label',
433                               '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id',
434                               '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos']
435                elif blk.get('_atom_site_moment_Fourier.atom_site_label'):
436                    MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label')
437                    MagFdict = dict(MagFloop.items())                           
438                    Mnames =  ['_atom_site_moment_fourier.atom_site_label',
439                               '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id',
440                               '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos']
441            self.Phase['Atoms'] = []
442            if magnetic:
443                self.MPhase['Atoms'] = []
444            G2AtomDict = {  '_atom_site_type_symbol' : 1,
445                            '_atom_site_label' : 0,
446                            '_atom_site_fract_x' : 3,
447                            '_atom_site_fract_y' : 4,
448                            '_atom_site_fract_z' : 5,
449                            '_atom_site_occupancy' : 6,
450                            '_atom_site_aniso_u_11' : 11,
451                            '_atom_site_aniso_u_22' : 12,
452                            '_atom_site_aniso_u_33' : 13,
453                            '_atom_site_aniso_u_12' : 14,
454                            '_atom_site_aniso_u_13' : 15,
455                            '_atom_site_aniso_u_23' : 16, }
456
457            ranIdlookup = {}
458            for aitem in atomloop:
459                atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
460                for val,key in zip(aitem,atomkeys):
461                    col = G2AtomDict.get(key,-1)
462                    if col >= 3:
463                        atomlist[col] = cif.get_number_with_esd(val)[0]
464                        if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
465                    elif col is not None and col != -1:
466                        atomlist[col] = val
467                    elif key in ('_atom_site_thermal_displace_type',
468                               '_atom_site_adp_type'):   #Iso or Aniso?
469                        if val.lower() == 'uani':
470                            atomlist[9] = 'A'
471                    elif key == '_atom_site_u_iso_or_equiv':
472                        uisoval = cif.get_number_with_esd(val)[0]
473                        if uisoval is not None: 
474                            atomlist[10] = uisoval
475                    elif key == '_atom_site_b_iso_or_equiv':
476                        uisoval = cif.get_number_with_esd(val)[0]
477                        if uisoval is not None: 
478                            atomlist[10] = uisoval/(8*np.pi**2)
479                if not atomlist[1] and atomlist[0]:
480                    typ = atomlist[0].rstrip('0123456789-+')
481                    if G2elem.CheckElement(typ):
482                        atomlist[1] = typ
483                    if not atomlist[1]:
484                        atomlist[1] = 'Xe'
485                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
486                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
487                    atomlist[9] = 'A'
488                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
489                        col = G2AtomDict.get(key)
490                        if col:
491                            atomlist[col] = cif.get_number_with_esd(val)[0]
492                if None in atomlist[11:17]:
493                    atomlist[9] = 'I'
494                atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
495                atomlist[1] = G2elem.FixValence(atomlist[1])
496                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
497                self.Phase['Atoms'].append(atomlist)
498                ranIdlookup[atomlist[0]] = atomlist[-1]
499                if atomlist[0] in atomlbllist:
500                    self.warnings += ' ERROR: repeated atom label: '+atomlist[0]
501                else:
502                    atomlbllist.append(atomlist[0])
503
504                if magnetic and atomlist[0] in magatomlabels:
505                    matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:]
506                    for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys):
507                        mcol = G2MagDict.get(mkey,-1)
508                        if mcol > 0:
509                            matomlist[mcol] = cif.get_number_with_esd(mval)[0]
510                    self.MPhase['Atoms'].append(matomlist)
511                if Super:
512                    Sfrac = np.zeros((4,2))
513                    Sadp = np.zeros((4,12))
514                    Spos = np.zeros((4,6))
515                    Smag = np.zeros((4,6))
516                    nim = -1
517                    waveType = 'Fourier'
518                    if  occCdict:
519                        for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']):
520                            if item == atomlist[0]:
521                                waveType = 'Crenel'
522                                val = occCdict['_atom_site_occ_special_func_crenel_c'][i]
523                                Sfrac[0][0] = cif.get_number_with_esd(val)[0]
524                                val = occCdict['_atom_site_occ_special_func_crenel_w'][i]
525                                Sfrac[0][1] = cif.get_number_with_esd(val)[0]
526                                nim = 1
527                   
528                    if nim >= 0:
529                        Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]]
530                    else:
531                        Sfrac = []
532                    nim = -1
533                    if displFdict:
534                        for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']):
535                            if item == atomlist[0]:
536                                waveType = 'Fourier'                               
537                                ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i])
538                                im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i])
539                                if im != nim:
540                                    nim = im
541                                val = displFdict['_atom_site_displace_fourier_param_sin'][i]
542                                Spos[im-1][ix] = cif.get_number_with_esd(val)[0]
543                                val = displFdict['_atom_site_displace_fourier_param_cos'][i]
544                                Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0]
545                    if nim >= 0:
546                        Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]]
547                    else:
548                        Spos = []
549                    nim = -1
550                    if UijFdict:
551                        nim = -1
552                        for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']):
553                            if item == atomlist[0]:
554                                ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i])
555                                im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i])
556                                if im != nim:
557                                    nim = im
558                                val = UijFdict['_atom_site_u_fourier_param_sin'][i]
559                                Sadp[im-1][ix] = cif.get_number_with_esd(val)[0]
560                                val = UijFdict['_atom_site_u_fourier_param_cos'][i]
561                                Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0]
562                    if nim >= 0:
563                        Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]]
564                    else:
565                        Sadp = []
566                    nim = -1
567                    if MagFdict:
568                        nim = -1
569                        for i,item in enumerate(MagFdict[Mnames[0]]):
570                            if item == atomlist[0]:
571                                ix = ['x','y','z'].index(MagFdict[Mnames[1]][i])
572                                im = int(MagFdict[Mnames[2]][i])
573                                if im != nim:
574                                    nim = im
575                                val = MagFdict[Mnames[3]][i]
576                                Smag[im-1][ix] = cif.get_number_with_esd(val)[0]
577                                val = MagFdict[Mnames[4]][i]
578                                Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0]
579                    if nim >= 0:
580                        Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]]
581                    else:
582                        Smag = []
583                    SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}}
584                    if magnetic and atomlist[0] in magatomlabels:
585                        matomlist.append(SSdict)
586                    atomlist.append(SSdict)
587            if len(atomlbllist) != len(self.Phase['Atoms']):
588                isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode'
589            for lbl in phasenamefields: # get a name for the phase
590                name = blk.get(lbl)
591                if name is None:
592                    continue
593                name = name.strip()
594                if name == '?' or name == '.':
595                    continue
596                else:
597                    break
598            else: # no name found, use block name for lack of a better choice
599                name = blknm
600            self.Phase['General']['Name'] = name.strip()
601            self.Phase['General']['Super'] = Super
602            self.Phase = copy.deepcopy(self.Phase)      #clean copy
603            if magnetic:
604                self.MPhase = copy.deepcopy(self.MPhase)    #clean copy
605                self.MPhase['General']['Type'] = 'magnetic'               
606                self.MPhase['General']['Name'] = name.strip()+' mag'
607                self.MPhase['General']['Super'] = Super
608                if Super:
609                    self.MPhase['General']['Modulated'] = True
610                    self.MPhase['General']['SuperVec'] = SuperVec
611                    self.MPhase['General']['SuperSg'] = SuperSg
612                    if self.MPhase['General']['SGData']['SGGray']:
613                        self.MPhase['General']['SuperSg'] += 's'
614                if 'mcif' not in filename:
615                    self.Phase = copy.deepcopy(self.MPhase)
616                    del self.MPhase
617            else:
618                self.MPhase = None
619            if Super:
620                if self.Phase['General']['SGData']['SGGray']:
621                    SGData = self.Phase['General']['SGData']
622                    SGData['SGGray'] = False
623                    ncen = len(SGData['SGCen'])
624                    SGData['SGCen'] = SGData['SGCen'][:ncen//2]
625                    self.Phase['General']['SGData'].update(SGData)
626                self.Phase['General']['Modulated'] = True
627                self.Phase['General']['SuperVec'] = SuperVec
628                self.Phase['General']['SuperSg'] = SuperSg
629                if not self.Phase['General']['SGData']['SGFixed']:
630                    self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1]
631            if self.ISODISTORT_test(blk):
632                if isodistort_warnings:
633                    self.warnings += isodistort_warnings
634                else:
635                    self.errors = "Error while processing ISODISTORT constraints"
636                    self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup)
637                    self.errors = ""
638            returnstat = True
639        return returnstat
640       
641    def ISODISTORT_test(self,blk):
642        '''Test if there is any ISODISTORT information in CIF
643
644        At present only _iso_displacivemode... and _iso_occupancymode... are
645        tested.
646        '''
647        for i in ('_iso_displacivemode_label',
648                  '_iso_occupancymode_label'):
649            if blk.get(i): return True
650        return False
651       
652    def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup):
653        '''Process ISODISTORT items to create constraints etc.
654        Constraints are generated from information extracted from
655        loops beginning with _iso_ and are placed into
656        self.Constraints, which contains a list of
657        :ref:`constraints tree items <Constraint_definitions_table>`
658        and one dict.
659        The dict contains help text for each generated ISODISTORT variable
660
661        At present only _iso_displacivemode... and _iso_occupancymode... are
662        processed. Not yet processed: _iso_magneticmode...,
663        _iso_rotationalmode... & _iso_strainmode...
664        '''
665        varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'}
666        'Maps ISODISTORT parm names to GSAS-II names'
667        # used for all types of modes
668        self.Constraints = []
669        explaination = {}
670        #----------------------------------------------------------------------
671        # read in the ISODISTORT displacement modes
672        #----------------------------------------------------------------------
673        if blk.get('_iso_displacivemode_label'):
674            modelist = []
675            shortmodelist = []
676            idlist = []
677            for id,lbl in zip(
678                blk.get('_iso_displacivemode_ID'),
679                blk.get('_iso_displacivemode_label')):
680                idlist.append(int(id))
681                modelist.append(lbl)
682                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
683            # just in case the items are not ordered increasing by id, sort them here
684            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
685            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
686            # read in the coordinate offset variables names and map them to G2 names/objects
687            coordVarLbl = []
688            G2varObj = []
689            idlist = []
690            error = False
691            for id,lbl in zip(
692                blk.get('_iso_deltacoordinate_ID'),
693                blk.get('_iso_deltacoordinate_label') ):
694                idlist.append(int(id))
695                coordVarLbl.append(lbl)
696                if '_' in lbl:
697                    albl = lbl[:lbl.rfind('_')]
698                    vlbl = lbl[lbl.rfind('_')+1:]
699                else:
700                    self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl
701                    error = True
702                    continue
703                if albl not in atomlbllist:
704                    self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl
705                    error = True
706                    continue
707                else:
708                    anum = atomlbllist.index(albl)
709                var = varLookup.get(vlbl)
710                if not var:
711                    self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl
712                    error = True
713                    continue
714                G2varObj.append(G2obj.G2VarObj(
715                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
716                    ))
717            if error:
718                raise Exception("Error decoding variable labels")
719            # just in case the items are not ordered increasing by id, sort them here
720            coordVarLbl = [i for i,j in sorted(zip(coordVarLbl,idlist),key=lambda k:k[1])]
721            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
722
723            if len(G2varObj) != len(modelist):
724                print ("non-square input")
725                raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate")
726
727            error = False
728            ParentCoordinates = {}
729            for lbl,exp in zip(
730                blk.get('_iso_coordinate_label'),
731                blk.get('_iso_coordinate_formula') ):
732                if '_' in lbl:
733                    albl = lbl[:lbl.rfind('_')]
734                    vlbl = lbl[lbl.rfind('_')+1:]
735                else:
736                    self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl
737                    error = True
738                    continue
739                if vlbl not in 'xyz' or len(vlbl) != 1:
740                    self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl
741                    error = True
742                    continue
743                i = 'xyz'.index(vlbl)
744                if not ParentCoordinates.get(albl):
745                    ParentCoordinates[albl] = [None,None,None]
746                if '+' in exp:
747                    val = exp.split('+')[0].strip()
748                    val = G2p3.FormulaEval(val)
749                    if val is None:
750                        self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl
751                        error = True
752                        continue
753                    ParentCoordinates[albl][i] = val
754                else:
755                    ParentCoordinates[albl][i] = G2p3.FormulaEval(exp)
756            if error:
757                print (self.warnings)
758                raise Exception("Error decoding variable labels")
759            # get mapping of modes to atomic coordinate displacements
760            displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj)))
761            for row,col,val in zip(
762                blk.get('_iso_displacivemodematrix_row'),
763                blk.get('_iso_displacivemodematrix_col'),
764                blk.get('_iso_displacivemodematrix_value'),):
765                displacivemodematrix[int(row)-1,int(col)-1] = float(val)
766            # Invert to get mapping of atom displacements to modes
767            displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix)
768            # create the constraints
769            modeVarList = []
770            for i,row in enumerate(displacivemodeInvmatrix):
771                constraint = []
772                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
773                    if k == 0: continue
774                    constraint.append([k,G2varObj[j]])
775                modeVar = G2obj.G2VarObj(
776                    (self.Phase['ranId'],None,shortmodelist[i],None))
777                modeVarList.append(modeVar)               
778                constraint += [modeVar,False,'f']
779                self.Constraints.append(constraint)
780            # normilization constants
781            normlist = []
782            idlist = []
783            for id,exp in zip(
784                blk.get('_iso_displacivemodenorm_ID'),
785                blk.get('_iso_displacivemodenorm_value'),
786                ):
787                idlist.append(int(id))
788                normlist.append(float(exp))
789            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
790            #----------------------------------------------------------------------
791            # save the ISODISTORT info for "mode analysis"
792            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
793            self.Phase['ISODISTORT'].update({
794                # coordinate items
795                'IsoVarList' : coordVarLbl,
796                'G2VarList' : G2varObj,
797                'ParentStructure' : ParentCoordinates,
798                # mode items
799                'IsoModeList' : modelist,
800                'G2ModeList' : modeVarList,
801                'NormList' : normlist,
802                # transform matrices
803                'Var2ModeMatrix' : displacivemodeInvmatrix,
804                'Mode2VarMatrix' : displacivemodematrix,
805                })
806            # make explaination dictionary
807            for mode,shortmode in zip(modelist,shortmodelist):
808                explaination[str(self.Phase['ranId'])+shortmode] = (
809                    "ISODISTORT full name "+str(mode))
810        #----------------------------------------------------------------------
811        # now read in the ISODISTORT occupancy modes
812        #----------------------------------------------------------------------
813        if blk.get('_iso_occupancymode_label'):
814            modelist = []
815            shortmodelist = []
816            idlist = []
817            for id,lbl in zip(
818                blk.get('_iso_occupancymode_ID'),
819                blk.get('_iso_occupancymode_label')):
820                idlist.append(int(id))
821                modelist.append(lbl)
822                ISODISTORT_shortLbl(lbl,shortmodelist) # shorten & make unique
823            # just in case the items are not ordered increasing by id, sort them here
824            modelist = [i for i,j in sorted(zip(modelist,idlist),key=lambda k:k[1])]
825            shortmodelist = [i for i,j in sorted(zip(shortmodelist,idlist),key=lambda k:k[1])]
826            # read in the coordinate offset variables names and map them to G2 names/objects
827            occVarLbl = []
828            G2varObj = []
829            idlist = []
830            error = False
831            for id,lbl in zip(
832                blk.get('_iso_deltaoccupancy_ID'),
833                blk.get('_iso_deltaoccupancy_label') ):
834                idlist.append(int(id))
835                occVarLbl.append(lbl)
836                if '_' in lbl:
837                    albl = lbl[:lbl.rfind('_')]
838                    vlbl = lbl[lbl.rfind('_')+1:]
839                else:
840                    self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl
841                    error = True
842                    continue
843                if albl not in atomlbllist:
844                    self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl
845                    error = True
846                    continue
847                else:
848                    anum = atomlbllist.index(albl)
849                var = varLookup.get(vlbl)
850                if not var:
851                    self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl
852                    error = True
853                    continue
854                G2varObj.append(G2obj.G2VarObj(
855                    (self.Phase['ranId'],None,var,ranIdlookup[albl])
856                    ))
857            if error:
858                raise Exception("Error decoding variable labels")
859            # just in case the items are not ordered increasing by id, sort them here
860            occVarLbl = [i for i,j in sorted(zip(occVarLbl,idlist),key=lambda k:k[1])]
861            G2varObj = [i for i,j in sorted(zip(G2varObj,idlist),key=lambda k:k[1])]
862
863            if len(G2varObj) != len(modelist):
864                print ("non-square input")
865                raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy")
866
867            error = False
868            ParentOcc = {}
869            for lbl,exp in zip(
870                blk.get('_iso_occupancy_label'),
871                blk.get('_iso_occupancy_formula') ):
872                if '_' in lbl:
873                    albl = lbl[:lbl.rfind('_')]
874                    vlbl = lbl[lbl.rfind('_')+1:]
875                else:
876                    self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl
877                    error = True
878                    continue
879                if vlbl != 'occ':
880                    self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl
881                    error = True
882                    continue
883                if '+' in exp:
884                    val = exp.split('+')[0].strip()
885                    val = G2p3.FormulaEval(val)
886                    if val is None:
887                        self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl
888                        error = True
889                        continue
890                    ParentOcc[albl] = val
891            if error:
892                raise Exception("Error decoding occupancy labels")
893            # get mapping of modes to atomic coordinate displacements
894            occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj)))
895            for row,col,val in zip(
896                blk.get('_iso_occupancymodematrix_row'),
897                blk.get('_iso_occupancymodematrix_col'),
898                blk.get('_iso_occupancymodematrix_value'),):
899                occupancymodematrix[int(row)-1,int(col)-1] = float(val)
900            # Invert to get mapping of atom displacements to modes
901            occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix)
902            # create the constraints
903            modeVarList = []
904            for i,row in enumerate(occupancymodeInvmatrix):
905                constraint = []
906                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
907                    if k == 0: continue
908                    constraint.append([k,G2varObj[j]])
909                modeVar = G2obj.G2VarObj(
910                    (self.Phase['ranId'],None,shortmodelist[i],None))
911                modeVarList.append(modeVar)               
912                constraint += [modeVar,False,'f']
913                self.Constraints.append(constraint)
914            # normilization constants
915            normlist = []
916            idlist = []
917            for id,exp in zip(
918                blk.get('_iso_occupancymodenorm_ID'),
919                blk.get('_iso_occupancymodenorm_value'),
920                ):
921                idlist.append(int(id))
922                normlist.append(float(exp))
923            normlist = [i for i,j in sorted(zip(normlist,idlist),key=lambda k:k[1])]
924            #----------------------------------------------------------------------
925            # save the ISODISTORT info for "mode analysis"
926            if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {}
927            self.Phase['ISODISTORT'].update({
928                # coordinate items
929                'OccVarList' : occVarLbl,
930                'G2OccVarList' : G2varObj,
931                'BaseOcc' : ParentOcc,
932                # mode items
933                'OccModeList' : modelist,
934                'G2OccModeList' : modeVarList,
935                'OccNormList' : normlist,
936                # transform matrices
937                'Var2OccMatrix' : occupancymodeInvmatrix,
938                'Occ2VarMatrix' : occupancymodematrix,
939                })
940            # make explaination dictionary
941            for mode,shortmode in zip(modelist,shortmodelist):
942                explaination[str(self.Phase['ranId'])+shortmode] = (
943                    "ISODISTORT full name "+str(mode))
944        #----------------------------------------------------------------------
945        # done with read
946        #----------------------------------------------------------------------
947        if explaination: self.Constraints.append(explaination)
948
949        if not debug: return
950
951        # debug: show displacive mode var to mode relations
952        if 'IsoVarList' in self.Phase['ISODISTORT']:
953            # coordinate items
954            #coordVarLbl = self.Phase['ISODISTORT']['IsoVarList']
955            G2varObj = self.Phase['ISODISTORT']['G2VarList']
956            #ParentCoordinates = self.Phase['ISODISTORT']['ParentStructure']
957            # mode items
958            modelist = self.Phase['ISODISTORT']['IsoModeList']
959            modeVarList = self.Phase['ISODISTORT']['G2ModeList']
960            normlist = self.Phase['ISODISTORT']['NormList']
961            # transform matrices
962            #displacivemodeInvmatrix = self.Phase['ISODISTORT']['Var2ModeMatrix']
963            #displacivemodematrix = self.Phase['ISODISTORT']['Mode2VarMatrix']
964       
965            print( 70*'=')
966            print('\nVar2ModeMatrix' ,'IsoVarList' )
967            for i,row in enumerate(displacivemodeInvmatrix):
968                l = ''
969                for j,(lbl,k) in enumerate(zip(coordVarLbl,row)):
970                    if k == 0: continue
971                    if l: l += ' + '
972                    #l += lbl+' * '+str(k)
973                    l += str(G2varObj[j])+' * '+str(k)
974                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
975
976            # Get the ISODISTORT offset values
977            coordVarDelta = {}
978            for lbl,val in zip(
979                blk.get('_iso_deltacoordinate_label'),
980                blk.get('_iso_deltacoordinate_value'),):
981                coordVarDelta[lbl] = float(val)
982            modeVarDelta = {}
983            for lbl,val in zip(
984                blk.get('_iso_displacivemode_label'),
985                blk.get('_iso_displacivemode_value'),):
986                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
987
988            print( 70*'=')
989            print('\nInverse relations using Var2ModeMatrix, NormList, IsoVarList')
990            # compute the mode values from the reported coordinate deltas
991            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
992                l = ''
993                for lbl,k in zip(coordVarLbl,row):
994                    if k == 0: continue
995                    if l: l += ' + '
996                    l += lbl+' * '+str(k)
997                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
998            print('\nCalculation checks\n')
999            for i,(row,n) in enumerate(zip(displacivemodeInvmatrix,normlist)):
1000                #l = ''
1001                sl = ''
1002                s = 0.
1003                for lbl,k in zip(coordVarLbl,row):
1004                    if k == 0: continue
1005                    #if l: l += ' + '
1006                    #l += lbl+' * '+str(k)
1007                    if sl: sl += ' + '
1008                    sl += str(coordVarDelta[lbl])+' * '+str(k)
1009                    s += coordVarDelta[lbl] * k
1010                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
1011                print(' ?= ',modeVarDelta[modelist[i]])
1012                print()
1013
1014            print( 70*'=')
1015            print('\nDirect relations using Mode2VarMatrix, NormList, IsoVarList')
1016            # compute the coordinate displacements from the reported mode values
1017            DeltaCoords = {}
1018            for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix):
1019                l = ''
1020                s = 0.0
1021                at,d = lbl.split('_')
1022                if at not in DeltaCoords:
1023                    DeltaCoords[at] = [0,0,0]
1024                for j,(k,n) in enumerate(zip(row,normlist)):
1025                    if k == 0: continue
1026                    if l: l += ' + '
1027                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
1028                    s += n * modeVarDelta[modelist[j]] * k
1029                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
1030                if d == 'dx':
1031                    DeltaCoords[at][0] = s
1032                elif d == 'dy':
1033                    DeltaCoords[at][1] = s
1034                elif d == 'dz':
1035                    DeltaCoords[at][2] = s
1036                else:
1037                    print('unexpected',d)
1038
1039            print( 70*'=')
1040            print('\nCoordinate checks')
1041            print('\nComputed')
1042            for at in sorted(DeltaCoords):
1043                s = at
1044                for i in range(3):
1045                    s += '  '
1046                    s += str(ParentCoordinates[at][i]+DeltaCoords[at][i])
1047                print(s)
1048
1049            # determine the coordinate delta values from deviations from the parent structure
1050            print('\nFrom CIF')
1051            for atmline in self.Phase['Atoms']:
1052                lbl = atmline[0]
1053                x,y,z = atmline[3:6]
1054                print( lbl,x,y,z)
1055
1056            print( 70*'=')
1057            print('\nGenerated constraints')
1058            for i in self.Constraints:
1059                if type(i) is dict:
1060                    print('\nconstraint help dict')
1061                    for key in i:
1062                        print('\t',key,':',i[key])
1063                elif i[-1] == 'f':
1064                    print('\n\t',i[-3],' =')
1065                    for m,j in i[:-3]:
1066                        print('\t\t+',m,' * ',j,'   ',repr(j))
1067                else:
1068                    print('  unexpected: ',repr(i))
1069            print("\nG2name ==> ISODISTORT full name",
1070                      ".Phase['ISODISTORT']['IsoModeList']",
1071                      ".Phase['ISODISTORT']['G2ModeList']")
1072            for mode,G2mode in zip(modelist,modeVarList):
1073                print("  ?::"+str(G2mode),' ==>', mode)
1074
1075        #======================================================================
1076        # debug: show occupancy mode var to mode relations
1077        if 'G2OccVarList' in self.Phase['ISODISTORT']:
1078            # coordinate items
1079            #occVarLbl = self.Phase['ISODISTORT']['OccVarList']
1080            G2varObj = self.Phase['ISODISTORT']['G2OccVarList']
1081            #ParentOcc = self.Phase['ISODISTORT']['BaseOcc']
1082            # mode items
1083            modelist = self.Phase['ISODISTORT']['OccModeList']
1084            modeVarList = self.Phase['ISODISTORT']['G2OccModeList']
1085            normlist = self.Phase['ISODISTORT']['OccNormList']
1086            # transform matrices
1087            #occupancymodeInvmatrix = self.Phase['ISODISTORT']['Var2OccMatrix']
1088            #occupancymodematrix = self.Phase['ISODISTORT']['Occ2VarMatrix']
1089           
1090            print( 70*'=')
1091            print('\nVar2OccMatrix' ,'OccVarList' )
1092            for i,row in enumerate(occupancymodeInvmatrix):
1093                l = ''
1094                for j,(lbl,k) in enumerate(zip(occVarLbl,row)):
1095                    if k == 0: continue
1096                    if l: l += ' + '
1097                    #l += lbl+' * '+str(k)
1098                    l += str(G2varObj[j])+' * '+str(k)
1099                print( str(i) + ': '+str(modeVarList[i])+' = '+l)
1100
1101            # Get the ISODISTORT offset values
1102            occVarDelta = {}
1103            for lbl,val in zip(
1104                blk.get('_iso_deltaoccupancy_label'),
1105                blk.get('_iso_deltaoccupancy_value'),):
1106                occVarDelta[lbl] = float(val)
1107            modeVarDelta = {}
1108            for lbl,val in zip(
1109                blk.get('_iso_occupancymode_label'),
1110                blk.get('_iso_occupancymode_value'),):
1111                modeVarDelta[lbl] = cif.get_number_with_esd(val)[0]
1112
1113            print( 70*'=')
1114            print('\nInverse relations using Var2OccModeMatrix, OccNormList, OccVarList')
1115            # compute the mode values from the reported coordinate deltas
1116            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
1117                l = ''
1118                for lbl,k in zip(occVarLbl,row):
1119                    if k == 0: continue
1120                    if l: l += ' + '
1121                    l += lbl+' * '+str(k)
1122                print('a'+str(i)+' = '+str(modeVarList[i])+' = ('+l+')/'+str(n))
1123            print('\nCalculation checks\n')
1124            for i,(row,n) in enumerate(zip(occupancymodeInvmatrix,normlist)):
1125                #l = ''
1126                sl = ''
1127                s = 0.
1128                for lbl,k in zip(occVarLbl,row):
1129                    if k == 0: continue
1130                    #if l: l += ' + '
1131                    #l += lbl+' * '+str(k)
1132                    if sl: sl += ' + '
1133                    sl += str(occVarDelta[lbl])+' * '+str(k)
1134                    s += occVarDelta[lbl] * k
1135                print(str(modeVarList[i]),'=','('+sl+') / ',n,'=',s/n)
1136                print(' ?= ',modeVarDelta[modelist[i]])
1137                print()
1138
1139            print( 70*'=')
1140            print('\nDirect relations using Occ2VarMatrix, OccNormList, OccVarList')
1141            # compute the coordinate displacements from the reported mode values
1142            Occ = {}
1143            for i,lbl,row in zip(range(len(occVarLbl)),occVarLbl,occupancymodematrix):
1144                l = ''
1145                s = 0.0
1146                for j,(k,n) in enumerate(zip(row,normlist)):
1147                    if k == 0: continue
1148                    if l: l += ' + '
1149                    l += str(n)+' * '+str(modeVarList[j])+' * '+str(k)
1150                    s += n * modeVarDelta[modelist[j]] * k
1151                print( lbl,'=',str(G2varObj[i]),'=',l,'=',s,'\n')
1152                j = lbl.split('_')[0] 
1153                Occ[j] = ParentOcc[j]+s
1154               
1155            # determine the coordinate delta values from deviations from the parent structure
1156            print('\nOccupancy from CIF vs computed')
1157            for atmline in self.Phase['Atoms']:
1158                lbl = atmline[0]
1159                if lbl in Occ: print( lbl,atmline[6],Occ[lbl])
1160
1161            print( 70*'=')
1162            print('\nGenerated constraints')
1163            for i in self.Constraints:
1164                if type(i) is dict:
1165                    print('\nconstraint help dict')
1166                    for key in i:
1167                        print('\t',key,':',i[key])
1168                elif i[-1] == 'f':
1169                    print('\n\t',i[-3],' =')
1170                    for m,j in i[:-3]:
1171                        print('\t\t+',m,' * ',j,'   ',repr(j))
1172                else:
1173                    print('  unexpected: ',repr(i))
1174            print("\nG2name ==> ISODISTORT full name",
1175                      ".Phase['ISODISTORT']['OccModeList']",
1176                      ".Phase['ISODISTORT']['G2OccModeList']")
1177            for mode,G2mode in zip(modelist,modeVarList):
1178                print("  ?::"+str(G2mode),' ==>', mode)
1179               
1180def ISODISTORT_shortLbl(lbl,shortmodelist):
1181    '''Shorten model labels and remove special characters
1182    '''
1183    regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)\[.*?\](.*)',lbl)
1184    # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)[...]BBBBB
1185    # where SSSSS is the parent spacegroup, [x,y,z] is a location, etc.
1186    if regexp:
1187        # this extracts the AAAAA and BBBBB parts of the string
1188        lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
1189    else:
1190        # try form SSSSS[x,y,z]AAAA(a,b,...)BBBBB
1191        regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl)
1192        if regexp:
1193            lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version
1194    lbl = lbl.replace('order','o')
1195    lbl = lbl.replace('+','_')
1196    lbl = lbl.replace('-','_')
1197    G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list
1198               
Note: See TracBrowser for help on using the repository browser.